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Abstract. Explicit computational formulas for coefficients of the periodic normal forms 
of the three most complex codim 2 bifurcations of limit cycles with dimension of the cen- 
ter manifold equal to 4 or to 5 in generic autonomous ODEs are derived. The resulting 
formulas are independent of the dimension of the phase space and involve solutions of 
certain boundary-value problems as well as multilinear functions from the Taylor expan- 
sion of the ODE right-hand side near the cycle. The formulas allow one to distinguish 
between the complicated bifurcation scenarios which can happen near these codim 2 
bifurcations of limit cycles, where 3-tori and 4-tori can be present. We apply our tech- 
niques to the study of a known laser model, a novel model from population biology, 
and one for mechanical vibrations; these models exhibit Limit Point-Neimark-Sacker, 
Period-Doubling-Neimark-Sacker and double Neimark-Sacker bifurcations. Lyapunov 
exponents are computed to numerically confirm the results of the normal form analysis, 
in particular with respect to the existence of stable invariant tori of various dimensions 
and chaos. 



1. Introduction 

Consider a smooth system of ODEs 

± = f(x,p), x£R n , (1.1) 

smoothly depending on a parameter vector p £ M m . Typically, the dynamics of such 
systems show qualitative transitions, i.e. bifurcations, upon variation of a parameter. 
It is hard to use simulations to characterize such transitions correctly and efficiently. 
Numerical continuation software such as auto [11] or matcont [10, 8, 9] may be used to 
track bifurcations from a stable equilibrium to a periodic oscillation by a Hopf bifurcation 
and even the appearance of (un)stable invariant tori with multi-frequency oscillations by a 
secondary Hopf, or Neimark-Sacker bifurcation. Bifurcations of these invariant tori Y m - 2 
into other tori or chaos, however, are out of reach of the standard numerical analysis. 

One possibility to study bifurcations of tori - if they are stable - is to compute Lyapunov 
exponents. The dimension of the torus for a given parameter value then equals the number 
of exponents equal to zero. Varying one parameter one can observe that exponents become 
zero and this indicates a bifurcation. The exact nature of the bifurcation is however 
obscured from this analysis and should be elucidated with additional means. Yet, in many 
cases, bifurcations of tori first emerge from codim 2 bifurcations of limit cycles. Specifically, 
these codim 2 bifurcations are points in the parameter plane where one Neimark-Sacker 
bifurcation curve intersects a Limit Point of cycles, a Period-Doubling or another Neimark- 
Sacker bifurcation curve. The intersections produce LPNS, PDNS, or NSNS bifurcations, 
respectively. This paper focuses on these bifurcations, occuring in generic systems (1.1) 
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when m > 2 and n is sufficiently large. The bifurcations are well understood theoretically 
with Poincare maps and the corresponding normal forms [3, 17, 15, 23, 27, 14, 33]. The 
results of the analyis of the normal form for these codim 2 bifurcations can be used to 
verify nondegeneracy conditions and classify the bifurcation structure. Hence, we need an 
algorithm for the numerical computation of the coefficients of each critical normal form 
to enable this analysis. 

There is a straightforward approach to obtain the critical normal forms of the codim 2 
bifurcations of the limit cycle. In the Poincare map, the limit cycle is a fixed point and 
one can use techniques developed for maps to obtain the critical normal form [27, 14]. 
However, in this case partial derivatives of the map up to order k, most often k = 3, 
sometimes k = 5, are needed. This may be done using software such as CAPD [1] or tides 
[2]. These packages can compute the solution and the derivatives of the solution with 
respect to the initial condition with arbitrary precision using Taylor series. Alternatively 
one could integrate the variational equations [ ] or use automatic differentiation [16, 26] 
to obtain the derivatives of the Poincare map. All these methods, however, have two 
drawbacks that make them less (time) efficient. First, these are shooting methods that 
are slower when the system is very sensitive to perturbations. Second, the full Poincare 
map is computed while only certain expressions are needed for the normalization. There 
is an alternative technique that is more suitable in the context of numerical continuation 
of periodic orbits using collocation as the whole periodic orbit is available. It uses periodic 
normalization [18, 19] and has been applied to codim 1 bifurcations of limit cycles and 
implemented in matcont [24]. Recently, we have extended this algorithm to codim 2 
bifurcations of limit cycles with center manifold dimension at most 3 [7] . Here we consider 
the three remaining cases, LPNS, PDNS, and NSNS, that are characterized by a center 
manifold of the critical cycle of dimension 4 or 5. These three cases always involve a - 
possibly unstable - two-dimensional torus T 2 . 

We have implemented our algorithm in the numerical continuation toolbox matcont 
which automatically invokes the algorithm whenever the corresponding bifurcation is de- 
tected. Hence, any user is able to use it and take advantage of the automated normal 
form analysis. Here we document precisely what our algorithm does. First, its aim is to 
compute coefficients of a periodic critical normal form. We present these normal forms in 
Section 2 using (contrary to [24, 7]) the original Iooss [18] representation. Remark that 
these normal forms are closely related to the normal forms for the Zero-Hopf and Hopf- 
Hopf bifurcations of equilibria. We discuss the correspondence and the interpretation of 
the bifurcation diagrams of the generic unfoldings for the LPNS, PDNS, and NSNS bifur- 
cations. Next, we present the formulas to compute the critical normal form coefficients in 
Section 3. Here we also comment on the implementation which is similar to [7]. Finally 
in Section 4, we consider several examples that involve tori bifurcations: a laser model, 
a model from population biology, and one for mechanical vibrations. In these models we 
find and analyze the three codim 2 bifurcations that we focus on. We compute the critical 
normal form coefficients using our algorithm to predict the bifurcation diagram near each 
of these codim 2 points. Next we corroborate the predictions using Lyapunov exponents. 
In fact, we argue that the classification from the critical normal form guides the correct 
interpretation of the Lyapunov exponents. 

2. Normal forms on the center manifold and their bifurcations 
Write (1.1) at the critical parameter values as 
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and suppose that there is a limit cycle T corresponding to a periodic solution uo(t) = 
uo(t + T), where T > is its (minimal) period. Expand F(uo(t) + v) into the Taylor series 

F(u (t) + v) = F(«o(t))+ 

A{t)v + -B{t- v, v) + -C(t; v, v, v)+ (2 .2) 
—D(t;v,v,v,v) + —E(t;v,v,v,v,v) +0(\\v\f), 

where A(t) = F u (uo(t)) and 

B(t;v 1 ,v 2 ) = F uu (u (t))[vi,v 2 ], C(t;v 1 ,v 2 ,v 3 ) = F uuu (u (t))[v 1 ,V2,v 3 ], 

etc. The matrix A and the multilinear forms B,C,D, and E are periodic in t with period 
T but this dependence will often not be indicated explicitly. 

Consider the initial- value problem for the fundamental matrix solution Y(t), namely, 

dY 

-£ = A{t)Y, Y(0) = I n , 

where I n is the nxn identity matrix. The eigenvalues of the monodromy matrix M = Y(T) 
are called (Floquet) multipliers of the limit cycle. The multipliers with = 1 are called 
critical. There is always a "trivial" critical multiplier fj, n = 1. We denote the total number 
of critical multipliers by n c and assume that the limit cycle is non- hyperbolic, i.e. n c > 1. 
In this case, there exists an invariant n c -dimensional critical center manifold W C (T) C W 1 
near r . 

2.1. Critical normal forms. It is well known [3, 23] that in generic two-parameter 
systems (1.1) only eleven codim 2 local bifurcations of limit cycles occur. To describe 
the normal forms of (2.1) on the critical center manifold IF c (r) for these codim 2 cases, 
we parameterize VF c (r) near T by (n c — 1) transverse coordinates and r G [0, kT) for 
k G {1,2,3,4}, depending on the bifurcation. The 8 cases where n c < 3 were treated in 
j7]. Based on [18] we show in Appendix A that the restriction of (2.1) to the corresponding 
critical center manifold M^ c (r) with n c = 4 or n c = 5 will take one of the following Iooss 
normal forms. 

2.1.1. LPNS. The Limit Point - Neimark-S acker bifurcation occurs when the trivial crit- 
ical multiplier fi n = 1 corresponds to a two-dimensional Jordan block and there are only 
two more critical simple multipliers fix^ = e ±ie with 9 ^ ^j-, for j = 1,2,3,4. The four- 
dimensional Iooss normal form at the LPNS bifurcation is derived in Appendix A. 1.1 and 
can be written as 

= 1 - 6 + "200^1 + "on |6| 2 + «30oCi + am£i + • • • , 

= a 2 oo£i + aon I6I 2 + a30o£i + amfi \&\ 2 + ■■■, (2.3) 

^ = iuj& + &no£i6 + &210C16 + &02i6 16I 2 + • • • , 

where r G [0, T], u> = ^1 is a real coordinate and £2 is a complex coordinate on IF c (r) 
transverse to T, a^^aijk G M, G C, and the dots denote the 0(||£ 4 ||)-terms, which 
are T-periodic in r. The equations (2.3) implicitly describe motions on the 4-dimensional 
invariant manifold lF c (r) with one cyclic coordinate r. 



^This manifold should not be confused with the (n c — l)-dimensional center manifold of the corresponding 
Poincare map. 
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2.1.2. PDNS. The Period- Doubling - Neimark-S acker bifurcation occurs when the trivial 
critical multiplier fi n = 1 is simple and there are only three more critical simple multipliers, 



namely —1 and /xi^ = e with 6 ^ -j-, for j = 1,2,3,4. The four-dimensional Iooss 



z ±ie with / ?f, 

normal form at the PDNS bifurcation is derived in Appendix A. 1.2 and can be written as 

( -^7 = 1 + "200^1 + "Oil |6| 2 + "40o£l + "022 |6| 4 + "21l£l |6| 2 + • • • > 

-j^- = O30o£i + oinCi |6| 2 + a 5 oo£i + 01226 I6| 4 + 03116* |6| 2 + • • • , 

^ = + 6210C16 + 60216 I6I 2 + &4io£i6 + &22i£i&2 I6I 2 + 60326 I6I 4 + • • • , 

(2.4) 

where r E [0, 2T], w = 0/T, £1 is a real coordinate and 6 is a complex coordinate on W C (T) 
transverse to T, ay^, ay& G R, foyfc € C, and the dots denote the 0(||£ 6 ||)-terms, which are 
2T-periodic in r. The equations (2.4) implicitly describe motions on the 4-dimensional 
invariant manifold W c (r) that is doubly covered by the selected coordinates. 

2.1.3. NSNS. The double Neimark-S acker bifurcation occurs when the trivial critical mul- 
tiplier [i n = 1 is simple and there are only four more critical simple multipliers = e ±i01 
and ^2,3 = e ±id2 with 6> lj2 / y , for j = 1, 2, 3, 4, 5, 6 and W 1 + j6 2 for l,j 6Z with Z+j < 4 
(see [14]). The five-dimensional periodic normal form at the NSNS bifurcation is derived 
in Appendix A. 1.3 and can be written as 

1 + "1100 |6| 2 + "0011 16| 2 



dt 



+"2200 |£l| 4 + "0022 |6| 4 + "1111 |£l| 2 I6| 2 + 



^ = iwi£i + a 2 iooCi |6| 2 + 010116 |6| 2 ^5) 
+032006 I6| 4 + 010226 |6| 4 + 021116 I6| 2 |6| 2 + • ■ ■ , 

iw 2 6 + 600216 |6l 2 + 611106 16I 2 



fir 



-600326 16I 4 + 622106 16I 4 + 611216 16I 2 16I 2 + 



where r E [0,T], u\ t 2 = Oip/T, 6 an d 6 are complex coordinates on T^ c (r) transverse to 
T, aijki E M.,aijki,bijki G C, and the dots denote the 0(||£ 6 ||)-terms, which are T-periodic 
in r. The equations (2.5) implicitly describe motions on a 5-dimensional manifold with 
one cyclic coordinate r. 

2.2. Generic unfoldings of the critical normal forms. Here we describe how the 
coefficients of the critical normal forms can be used to predict bifurcations of the phase 
portraits near the critical limit cycles for nearby parameter values. We introduce certain 
quantities - computable in terms of these coeffcients - that are reported in the matcont 
output and used to distinguish between various bifurcation scenarios in examples in Section 
4. 

In generic two-parameter systems (1.1) the considered bifurcations occur at isolated 
parameter values. By translating the origin of the parameter plane to one of such points, 
we can consider an unfolding of the corresponding bifurcation and study its canonical 
local bifurcation diagram for nearby parameter values. It is well known that the critical 
center manifold H^ c (r) can be smoothly continued w.r.t. p in a neighborhood of the 
bifurcation point, so that the restriction of (1.1) to this manifold can be studied. Choosing 
appropriate coordinates (6 t) on this parameter-dependent invariant manifold, one can 
transform the restricted system into a parameter-dependent normal form in which 4* has 
a r-independent principle part and higher-order terms which are fcT-periodic in r with 
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k = 1 for LPNS and NSNS and k = 2 for PDNS. Below we describe bifurcations of these 
principle parts, i.e., the truncated parameter-dependent autonomous normal forms. Since 
the dynamics is determined by the ^-equations, we first focus on their bifurcations and 
then interpret appearing bifurcation diagrams for the original system (1.1). The new 
unfolding parameters will be denoted by /3a). 

2.2.1. LPNS. Generically, a two-parameter unfolding of (1.1) near this bifurcation re- 
stricted to the center manifold is smoothly orbitally equivalent (with possible time rever- 
sal) to a system in which the equations for the transverse coordinates have the form 

§- = /3i+£ 2 + S |C| 2 + 0(||(£,C,C)H 4 ), 

Z (2-6) 

^ = ( / 3 2 + ^ 1 )c + ^ + ^)ec + e 2 c + o(ii(e,c,c)ii 4 ), 

where the O-terms are still T-periodic in r. This system is similar to the normal form for 
the Zero-Hopf bifurcation of equilibria (cf. Theorem 8.6 on page 338 in [23]). In Figure 
1 the four possible bifurcation diagrams of the amplitude system for (2.6) without the 
O-terms, 

p = p({3 2 + et; + e), [ ] 

are reported depending on the sign of the normal form coefficients s and 9 [23]. Here and 
in what follows a dot means the derivative w.r.t. r. 

Let us now discuss the interpretation of the phase portraits in the (£, /?)-plane of the 
truncated amplitude system in the context of the bifurcating limit cycle. The fixed points 
or limit cycles have additional dimensions from the phases of the periodic orbit itself plus 
the phases ignored in the reduction to the amplitude system. We note that in the amplitude 
system the vertical direction always corresponds to a Neimark-S acker bifurcation, but 
that the horizontal component of the phase space has a different meaning. For LPNS, 
equilibria on the horizontal axis correspond to limit cycles. Equilibria off the horizontal 
axis correspond to invariant 2D tori T 2 and the periodic orbit which exists if s9 < 
corresponds to an invariant 3D torus T 3 . 

The critical values of s and 9 can be expressed in terms of the coeffcients of (2.3) as 

s = sign (a 2 ooooii), 9 = . 

0200 

These values determine the bifurcation scenario. For s9 < 0, a 3-torus appears in the 
unfolding via a Neimark-S acker bifurcation. The stability of this torus is determined 
by the third order terms in (2.3). Indeed, the sign of the corresponding first Lyapunov 
coefficient for the Hopf bifurcation in (2.7) is opposite to that of 9 but the 'time' in (2.6) 
is rescaled with factor 

^ B /, . , ( 3£(&02l) 3a 30 o , am \ ^0210200 

£j = U[ &210 + 0110 7, r 

\ V a 011 ^«200 ^«011 / «011 

(see page 337 in [23]). If E ■ l\ < 0, an unstable 3-torus appears, if E • l\ > 0, the 3-torus 
is stable. The output given by matcont is (s, 9, E) 2 . 

Note that Figure 1 presents bifurcations of the truncated system (2.6) that only approxi- 
mates the full normalized unfolding. In particular, the orbit structure on the invariant tori 
can differ from that for the approximating system due to phase locking. Moreover, the de- 
struction of T 3 via a heteroclinic bifurcation in case (c) of Figure 1 becomes a complicated 
sequence of global bifurcations involving stable and unstable invariant sets of cycles and 



Remark that E = NaN is reported when terms up to only the second order are computed. 
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(c) s = 1,0 < (d) s = -1, > 



Figure 1. Bifurcation diagrams of the truncated amplitude system (2.7) 
for the LPNS bifurcation. 

tori. All these bifurcations, however, occur in the exponentially-small parameter wedge 
near the heteroclinic bifurcation curve P. For detailed discussions of the effects of the 
truncation, also in the two other cases, we refer to [27, 33] and references therein. 

2.2.2. PDNS. Generically, a two-parameter unfolding of (1.1) near this bifurcation re- 
stricted to the center manifold is smoothly orbitally equivalent to a system in which the 
equations for the transverse coordinates have the form 

( v l = n lVl + p uV f + P 12 vi \v 2 \ 2 + Sivi \v 2 \ 4 + 0(11^^2, v 2 )\\ 6 ), 
I i>2 = (fJ-2 + ioj 2 )v 2 + P2iv\v 2 + P 22 v 2 \v 2 \ 2 + S 2 vjv 2 + iR 2 v 2 \v 2 \ 4 + 0(\\(v 1 ,v 2 ,v 2 )\\ G ), 

(2.8) 

where the O-terms are still T-periodic in r. This system is similar to one of the normal 
forms for the Hopf-Hopf bifurcations of equilibria (cf. Lemma 8.14 on page 354 in [23]). 
The amplitude system for (2.8) without the O-terms is 



\r 2 = ri(n2+P2ir 2 +P22r2 + s 2 rf), 



(2.9) 
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where 

PU = Pll, Pl2 = Pl2, P21 = ^(P 2 l), p 2 2 = 5?(P 22 ), Si = St, S 2 = $t(S 2 ). 

The values of pjk and Sj, for j, k = 1,2, and the quantities 

9 = P -^,5 = ^,Q = ^,A = ^- 

P22 Pll P A 22 Pll 

indicate in which bifurcation scenario we are (see Section 8.6.2 in [23]). 

In the "simple" case where P11P22 > 0, there are five topologically different bifurcation 
diagrams of the truncated amplitude system (2.9), corresponding to the following cases: 
I. > 0, 5 > 0, 05 > 1 
II. 9 > 0, 5 > 0, 95 < 1 

III. 9 > 0, 5 < 

IV. 9 < 0, 5 < 0, 95 < 1 
V. 9 < 0, 5 < 0, 95 > 1 

If 5 > 9, reverse the role of 9 and 5. Each case corresponds with a region in the (9, 5)- 
plane, see Figure 2 (a). The parametric portraits belonging to the different regions can be 
seen in Figure 3 (a), with corresponding phase portraits in the (7*1, r2)-plane in Figure 3 
(b). The phase portraits are only shown for the case when pn < and P22 < 0. The case 
pu > and P22 > can be reduced to the considered one by reversing time. 

In the "difficult" case where P11P22 < however, there are six essentially different 
bifurcation diagrams: 

I. 9 > 1,5 > 1 
II. 6 > 1,5 < 1,95 > 1 

III. 9 > 0,5 > 0,95 < 1 

IV. 9 > 0, 5 < 

V. < 0, 5 < 0, 95 < 1 
VI. 9 < 0, 5 < 0, 95 > 1 
The regions in the (9, £)-plane are shown in Figure 2 (b). The related parametric portraits 
and phase portraits of (2.9) are given in Figure 4. Only the case pn > and P22 < is 
presented, to which the opposite one can be easily reduced. 

We note that Section 8.6.2 in [23] for the "difficult" case contains a few errors in the 
figures and in the asymptotic expression for the heteroclinic bifurcation curve 3 . Therefore, 
for completeness, we provide the figures and correct asymptotics in Appendix B. 

The critical values of Pjk and Sj can be expressed in terms of the coeffcients of (2.4) as 

P11 = 02100, P12 = a 1011 , K(P 21 ) = ^(6 mo ), 3?(P 22 ) = K(6 00 2i), 



and 



, '^(&H2i) «, ^(60032) 03200^(60021! 

<5l — O1022 + OlOll T^Tr \" — l \ 



K(5 2 ) = 3?(&22io) + K(6 



K(6mo) K(&002l) «2100^(^lllo) 
^2111 „ ^3200 02100^(^0032) 



1110 J 



aion a 2 ioo oionK(6oo2i) 
(see page 356 in [23]). 

The fifth-order terms in (2.4) determine the stability of the tori in the "difficult" cases. 
In fact, the sign of the first Lyapunov coefficient for the Neimark-S acker bifurcation is 
given by 

sign h = -sign {9(0(9 - 1) A + 5(5 - 1)6)) . (2.10) 



Unfortunately, there is also a minor misprint in our earlier "correction" for the heteroclinic curve given 
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(a) (b) 



Figure 2. (a) the five subregions in the 5)-plane in the "simple" case; 
(b) the six subregions in the (9, 5)-plane in the "difficult" case. 




(a) (b) 



Figure 3. Bifurcation diagrams of the amplitude system (2.9) for the 
PDNS and NSNS bifurcations: (a) parametric portraits in the "simple" 
case; (b) phase portraits in the "simple" case. 

The output of matcont is (pu,P22, 6, 6, sign ii) 4 . 



Remark that sign l\ = NaN is reported when terms up to only the third order are computed. 
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(a) (b) 



Figure 4. Bifurcation diagrams of the amplitude system (2.9) for the 
PDNS and NSNS bifurcations: (a) parametric portraits in the "difficult" 
case; (b) phase portraits in the "difficult" case. 

For PDNS we have an interpretation analogous to LPNS, but the invariant sets may 
be "doubled". The origin always corresponds the original limit cycle. Other fixed points 
on the horizontal axis represent the period-doubled limit cycles, while a fixed point on 
the vertical axis corresponds to a T 2 . Fixed points off the coordinate axes correspond to 
doubled tori T 2 and periodic orbits correspond to T 3 . As in the LPNS case, Figures 3 and 
4 present bifurcations of the truncated amplitude system that only approximates the full 
normalized unfolding. In particular, one has to be carefull with 'torus doubling', which is 
in fact a complicated quasiperiodic bifurcation [29, 32]. 

2.2.3. NSNS. Generically, a two-parameter unfolding of (1.1) near this bifurcation re- 
stricted to the center manifold is smoothly orbitally equivalent to a system in which the 
equations for the transverse coordinates have the form 

\vi = (hi + iui)v\ + Px\V\ \v\\ 2 + P\2V\ \v 2 \ 2 + iR\v\ \v\\ A + S\V\ \v 2 \ A + 0(\\(v, v)\\ 6 ), 

\i>2 = (/-i2 + ^2)^2 + P21V2 \vi\ 2 + P22V2 I ^2 1 2 + S2V2 |vi| 4 + 1R2V2 1 4 + 0(||(f,u)|| 6 ), 

(2.11) 

where the O-terms are T-periodic in r. Neglecting this periodicity, system (2.11) is the 
normal form for the Hopf-Hopf bifurcation of equilibria (cf. Lemma 8.14 on page 354 in 
[23]). 
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The truncated amplitude system for (2.11) is the same as (2.9), where now 
p u = »(Pn) = &(a 2 ioo), P12 = K(Pi 2 ) = »(oioii), 

P21 = K(P 2 l) = K(6nio), P22 = K(P 22 ) = 3?(6o02l), 



and 




The output of matcont is (pn,p 22 , #, <5, sign Zi) 5 . 

Although the phase portraits of the truncated amplitude system are the same as for 
PDNS, their interpretation is slightly different, since they 'live' in the (|t»i|, |v 2 |)-plane. 
Here, on both axes the fixed points correspond to invariant 2D tori T 2 for the original 
system. Fixed points off the coordinate axes and limit cycles correspond to T 3 and T 4 , 
respectively. The usual remark on the approximate nature of the bifurcation diagrams 
applies here as well. 



As was mentioned in the previous section, the stability of the extra torus appearing in 
the "difficult" cases is determined by third order terms for the LPNS bifurcation and fifth 
order terms for the PDNS and NSNS bifurcations. In the "simple" cases, second order 
derivatives are sufficient to determine the behaviour in the LPNS bifurcations and third 
order derivatives are sufficient in the PDNS and NSNS bifurcations. Therefore, we restrict 
our computations in this section to second order terms in the LPNS case and up to and 
including third order terms in the PDNS and NSNS cases. The expressions of the third 
order coefficients for LPNS and fourth and fifth order coefficients for PDNS and NSNS 
are given in Appendix C. Remark that for efficiency reasons these higher order coefficients 
are not computed in matcont, unless explicitly requested by the user. 

3.1. LPNS. The four-dimensional critical center manifold VF c (r) at the LPNS bifurcation 
can be parametrized locally by (£i,£ 2 ,t) elxCx [0, T] as 



3. Computation of critical coefficients 



u = u (t) + £iVi(t) + &V2(t) +C2^ 2 (r) +H(£i,£ 2 ,t), 
where H satisfies if(£i,£ 2 ,T) = H(^i,^ 2 ,0) and has the Taylor expansion 



(3.1) 



„ , . . ... ,„ o.l .K. 



(3.2) 




(3.3) 



and 



v 2 - A(t)v 2 + iujv 2 = 0, tG[0,T], 
V2(T)-v 2 (0) = 0, 
, Jo (V2,v 2 )dr - 1 = 0. 



(3.4) 



'Remark that sign l\ = NaN is reported when terms up to only the third order are computed. 
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The functions V\ and v 2 exist because of Lemma 2 of [18]. The functions hijk will 
be found by solving appropriate BVPs, assuming that (2.1) restricted to W c (r) has the 
normal form (2.3). 

The coefficients of the normal form arise from the solvability conditions for the BVPs as 
integrals of scalar products over the interval [0, T\. Specifically, those scalar products in- 
volve among other things the quadratic and cubic terms of (2.2) near the periodic solution 
uq, the generalized eigenfunction v% and eigenfunction v 2 , and the adjoint eigenfunctions 
<p*, v* and t>2 as solutions of the problems 

<p* + A T (T)<p* 
p*(T)-p*(0) 



Jo ' (<P*,V!)dl 



1 



0, r e [0,T] 
0, 
0. 



and 



f vf + A t {t)vI + if* 
vt(T)-v*M 
Jo(vl,vi)dT 

)l + A t (t)v% + iuv% 
v*(T)-v* 2 (0) 



0, r 6 [0, T], 

0, 

0, 

0, re [0,T], 

o, 

0. 



Jo (v 2 ,v 2 )dT-l 
In what follows we will make use of the orthogonality condition 

(-T 

{<p*,F(u o ))dr = 0, 



and the normalization condition 



(«t,F(«o))dr = l J 



(3.5) 



(3.6) 



(3.7) 



(3.8) 



(3.9) 



which can be easily obtained from (3.3), (3.5) and (3.6). 

To derive the normal form coefficients we write down the homological equation and 
compare term by term. We therefore substitute (3.1) into (2.1), using (2.2), (2.3) and 
(3.2). By collecting the constant and linear terms we get the identities 

u = F(u ), v± - F(uo) = A(t)vi, v 2 + iujv 2 = A(t)v 2 , 

and the complex conjugate of the last equation. 

By collecting the £^-terms we find an equation for h 2 $o 

h 2 oo ~ A{T)h 2m = B(t; vi,v{) - 2a 200 v 1 - 2a 2 o u + 2wi, (3.10) 

to be solved in the space of functions satisfying h 2 oo(T) = h 2 oo(0). In this space, the 
differential operator Jp— A(r) is singular and its null-space is spanned by uq. The Fredholm 
solvability condition 

rT 

(ip*,B(T;v 1 ,v 1 ) - 2a 200 vi - 2a 200 u + 2vi) dr = 

allows one to calculate the coefficient 0200 i n (2-3) due to the required normalization in 
(3.5), i.e. 



0-200 



1 r 



B(t; vi,vi) + 2A(t)vi) dr, 



(3.11) 



taking (3.3) and (3.8) into account. With 0200 defined in this way, let /1200 be a solution 
of (3.10) in the space of functions satisfying /i20o(0) = h 2 oo(T). Notice that if h 2 oo is a 
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solution of (3.10), then also h 2 oo + £iF(uo) satisfies (3.10), since F(uo) is in the kernel of 
the operator 4- — A(r). In order to obtain a unique solution (without a component along 
the null eigenspace) we impose the following orthogonality condition which determines the 
value of E\ 

cT 

(v*,h 2 oo) dr = 0, 

'o 

since (3.9) holds. Thus h 2 oo is the unique solution of the BVP 

^200 - A(r)h 2 oo - B(t; vi, vt) - 2A(t)v\ + 2o 20 ovi + 2a 2 oo'"o - 2tt = 0, r G [0, T], 

h 20 o(T) - 7»2oo(0) = 0, 
fo{v{,h 2 oo)dT = 0. 



f 

Jo 



By collecting the £|-terms (or £|-terms) we find an equation for /1020 
^020 - A(T)h 020 + 2iujh 02 o = B(t; v 2 ,v 2 ), 



(3.12) 



(or its complex conjugate). This equation has a unique solution /1020 satisfying ho 2 o(T) = 
^020(0), since due to the spectral assumptions e 2tulT is not a multiplier of the critical cycle. 
Thus, /1020 can be found by solving 



(3.13) 



^020 - A(T)h 020 + 2iw/i 2o - B(t; v 2 ,v 2 ) = 0, r G [0, T], 
h 020 {T) - h 020 {0) = 0. 

By collecting the ^1^2-terms we obtain an equation for huo 

huo - A(r)hno + iuhno = B ( T i n,v 2 ) - b lw v 2 + v 2 + iuv 2 , 

to be solved in the space of functions satisfying hno(T) = /iiro(0). In this space, the 
differential operator ^ — A(r) + iu) is singular, since e lU)T is a critical multiplier. So we 
can impose the following Fredholm solvability condition 

T 

(v 2 ,B(t;vi,v 2 ) - b no v 2 + v 2 + iujv 2 ) dr = 0, 

which due to the normalization condition in (3.7) determines the value of the normal form 
coefficient buo, yielding 



'iio= / (v 2 ,B(t;vi,v 2 ) + A(t)v 2 ) dr. 
Jo 



(3.14) 



The nullspace belonging to the operator — A(r) + iu is one-dimensional and spanned 
by v 2 . To determine huo uniquely, we need to impose an orthogonality condition with a 
vector whose inproduct with v 2 is non-zero. v 2 can be choosen because of the normalisation 
condition in (3.7). In fact, huo only appears in the normal form coefficient 6210 ( see 
Appendix C.l), and a different normalization of huo does not influence the value of that 
normal form coefficient. Therefore, we obtain huo as the unique solution of the BVP 

huo- A(r)huo+iuhuo- B(T;vi,v 2 ) + buov 2 - A(t)v 2 = 0, te[0,T], 

h lw (T) - h no (0) = 0, (3.15) 
Jo ( u 2 A10} dr = 0. 

By collecting the |^2| 2 -terms we obtain a singular equation for /ton, namely 



h u ~ A{r)hou = B(t; v 2 ,v 2 ) - aoil^l - "oil^o, 
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to be solved in the space of functions satisfying hon(T) = /ioii(0). The non-trivial kernel 
of the operator Jp— A(t) is spanned by uq. So, the following Fredholm solvability condition 
is involved 

rT 

(ip* ,B(t;v 2 ,v 2 ) - a on vi - aoniio) dr = 0, 

'o 

which gives us the expression for the normal form coefficient aon, i.e. 



((P*,B(t;v 2 ,v 2 )) dr. 



(3.16) 



We impose the orthogonality condition with the adjoint generalized eigenfunction to 
obtain /ion as the unique solution of 



' h on - A(T)h 0U - B(t;v 2 ,v 2 ) + aonV! + a 01 iu = 0, r G [0, T] 

h 011 (T) - h on (0) = 0, 
!q{v\Mii) dr = 0. 



(3.17) 



We remark that the values of a 2 oo and aon are n °t determined by the homological 
equation. We therefore put them equal to zero. 

Third order coefficients are only needed to determine the stability of the torus, if this 
torus exists. For completeness, we have listed these terms in Appendix C. 

3.2. PDNS. The four-dimensional critical center manifold VF c (r) at the PDNS bifurca- 
tion can be parametrized locally by (£i, £2; r) G R x C x [0, 2T] as 



u = u (t) +£l«i(t) +&v 2 (t) +i 2 v 2 {r) +H(£i,&,t), 
where H satisfies H(^i,^ 2 ,2T) = H(^i,^ 2 ,0) and has the Taylor expansion 



E 



1 



i\j\k\ 



hijk(T)a&i+o(u\\ 



(3.18) 



(3.19) 



2<i+i+fc<5 

while the eigenfunctions v i and v 2 are defined by 

v\ — A(t)v i = 
ui(T) + ui(0) = 
Jo ( v l,Vl)dT -1 = 

with v\{t + T) = — vi(t) for r € [0, T] and 

i) 2 - A(t)v 2 + iuiv 2 
v 2 (T) - v 2 (0) 
Jo (v 2 ,v 2 )dr -1 



0, re [0,T], 

0, 

0. 



0, t G [0, T], 

o, 

0. 



(3.20) 



(3.21) 



The functions v\ and v 2 exist because of Lemma 5 of [18]. The functions can 
be found by solving appropriate BVPs, assuming that (2.1) restricted to VF c (r) has the 
normal form (2.4). Moreover, n(r, £i, £ 2 , £ 2 ) = u(t + T, — £i, ^2, C2) so that 



hijkir) = (-i)% jk (r + T), 



(3.22) 



for r G [0, T]. Therefore, we can restrict our computations to the interval [0,T] instead of 
[0,2T]. 

The coefficients of the normal form arise from the solvability conditions for the BVPs as 
integrals of scalar products over the interval [0, T]. Specifically, those scalar products in- 
volve among other things the quadratic and cubic terms of (2.2) near the periodic solution 
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1*0) vi, v 2 , and the adjoint eigenfunctions 92*, v\ and v\ as solutions of the problems 

<p* + A T (r)<p* = 0, r€ [0,T], 

tp*(T) — ip*(0) = 0, (3.23) 
J T (<p*,F(u ))dT-l = 0, 

vI + A t {t)vI = 0, r G [0,T], 

uj(r)+«i(°) = °> (3-24) 
/ T (^,^}dT-l = 0, 

and 

?>| + A t (t>2* + iuvl = 0, re[0,T], 

v|(T)-^(0) = 0, (3.25) 

J T (^,« 2 )dr-l = 0. 

By collecting the constant and linear terms we get the identities 

iio = F(u ), v\ = A(t)vi, i> 2 + 10JV2 = A(t)v 2 , 

and the complex conjugate of the last equation, which merely reflect the definition of uq 
and (3.20), (3.21). 

By collecting the ^-terms we find an equation for H 2 qq 

h 2 oo - A(T)h 2 oo = B(t; vi,vi) - 2a 2 oo^o, (3.26) 

to be solved in the space of functions satisfying h 2 oo(T) = ^200 (0)- In this space, the 
A. 

dr' 

solvability condition 



differential operator 4-—A(j) is singular and its null-space is spanned by uq. The Fredholm 



/ (<p*, B(t; vx, vi) - 2a 2 oo«o) dr = 
Jo 



gives us the possibility to calculate 0200 in (2.3) by the required normalization in (3.23), 
i.e. 

1 f T 

"200 = 7:/ {w*,B{t;vi,V])) dr. (3.27) 



2 jo 

As before, /i2oo is determined up to the addition of a multiple of ?io, since /1200 + £\F{uq) 
is a solution of (3.26) for every value of £\. We fix the value of /1200 by demanding the 
orthogonality with the adjoint eigenfunction corresponding with multiplier 1, i.e. 

T 

h 20 o) dr = 0. 



We obtain /1200 then as the unique solution of the BVP 

h 2 oo - A(T)h 20 o - B(t;vi,vi) + 2a 2 oou = 0, rG[0,T], 

h 200 (T) - h 200 (0) = 0, (3.28) 
Jo (vA^oo) dr = 0. 
By collecting the £|-terms (or ^-terms) we find the differential equation for /1020 
ho20 ~ A(T)h 020 + 2iujh 020 = B(t; v 2 ,v 2 ), 

or its complex conjugate. Since e 2lLoT is not a critical multiplier, no Fredholm solvability 
condition has to be satisified. /2.020 can thus simply be found by solving 

ho20 ~ A(T)h 020 + 2iuih 020 - B(r; v 2 ,v 2 ) = 0, r <E [0, T], ,„ , 

^020(^-^020(0) = 0. [6 - zy> 
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The equation found by comparing the £i£2-terms is given by 
huo - A(T)h 11Q + iuhno = B(r; v 1 ,v 2 ). 

From (3.22) it follows that huo is anti-periodic. Now, since —e tU)T is not a multiplier of the 
critical cycle, no solvability condition has to be satisfied. Therefore, we can immediately 
obtain huo from 

h lw - A(T)h lw + iujhno - B(t;vi,v 2 ) = 0, rG[0,T], , . 

h 110 (T) + h uo (0) = 0, {6 - 6U) 

The |^2| 2 -terms lead to a singular equation for /ion, namely 

A ii - A(r)h u = B(t; v 2 ,v 2 ) - a ntto, 

to be solved in the space of T-periodic functions. The non-trivial kernel of the operator 
4p — A(t) is spanned by uq. So, the Fredholm solvability condition with the corresponding 
T-periodic adjoint eigenfunction is involved, i.e. 

T 

((p*,B(r; v 2 ,v 2 ) - a on uo) dr = 0, 



from which the expression for the normal form coefficient aon can be derived 

"011= / (<p*,B(t;v 2 ,v 2 )) dr. (3.31) 
J 

Now, we still need to uniquely determine the multiple of F(uo) which can be added to 
the function /ion, and will therefore impose the orthogonality condition with ip* to obtain 
/ion as the unique solution of 

h on - A(T)h on - B(t;v 2 ,v 2 ) + aonuo = 0, rG[0,T], 

h n(T)- /lon(O) = 0, (3.32) 
fa{<P*,hau)dT = 0. 

We have now examined all order two terms, and continue with the order three terms. 

Collecting the £f-terms gives an equation for /1300 and will give us the possibility to 
compute the normal form coefficient 0300 in (2.4). The differential equation 

A300 - A(T)h 300 = C(r; v\, v\, V\) + 3B(t; vi,h 200 ) - 6a 2 oo*i - 6a 30 o«i 

has to be solved in the space of functions satisfying /i3oo(T) = — ^300 (0)- The non-trivial 
anti-periodic kernel of the operator J^— A(r) is spanned by v\. So, the Fredholm solvability 
condition with the anti-periodic adjoint eigenfunction v* is involved, i.e. 

/ («*,C(r; vi,vi,vi) + 3B(t; vi, h 200 ) -60200^1 -6a 3 oo^i) dr = 
J 

and thus 



T 



0300 = g / (v*,C(T;v 1 ,vi,vi) + SB(r;vi,h2oo)-Qoc2ooA(r)vi) dr, 







(3.33) 



due to the normalization condition from (3.24). The usual orthogonality condition with 
the adjoint eigenfunction u* is imposed to obtain hsoo as the unique solution of 

A300 - A(r)h 30Q - C(r; vi,vi,vi) - 3B{r;vi,h 2 oo) 

+6a 200 A(T)v 1 + 6a 300 v 1 = 0, rG[0,T], 

h 300 (T) + h 300 (0) = 0, [6 - 64) 

f {v*,h 300 ) dr = 0. 
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The £2 ( or Cf)"terms from the homological equation give an equation for /1030 
^030 - A( T ) h 030 + 3iw/i 30 = C(t;v 2 ,v 2 ,v 2 ) + 3B(t;v 2 , h 020 ), 

or its complex conjugate. This equation has a unique solution ^030 satisfying ho3o(T) = 
^030 (0), since due to the spectral assumptions e 3tujT is not a multiplier of the critical cycle. 
Thus, /1030 can be found by solving 

h 30 ~ M T ) h 030 + 3iL}h 030 - C(t;v 2 ,v 2 ,v 2 ) - 3B(T-,v 2 ,h 020 ) = 0, rG[0,T], 

h 030 {T) - h 030 {0) = 0. 

(3.35) 

By collecting the ^^-terms we find an equation for h 2 ±o 

h 2W - A(T)h 2W + iujh 210 = C(r; v\, v 1} v 2 ) + B(t; v 2 , h 2 oo) + 2B(t; vt, hu ) , g 3g x 

- 20200^2 - 26 2 io«2 - 2iu)a 20 ov 2 , 

to be solved in the space of T-periodic functions. The non-trivial kernel of the operator 
Jp — A(t) + iu is spanned by the complex eigenfunction v 2 . So, the following Fredholm 
solvability condition has to be imposed 

T 

(vI,C{t;vx,vx,v 2 ) + B(T;v 2 ,h 200 ) + 2B(t;vi, hu ) 



-20200^2 - 26210^2 - 2iuja 200 v 2 ) dr = 0. 
From this, the expression for the normal form coefficient 6210 can be derived, namely 



1 f T 

&210 = o / \v 2 ,C(t;vi,vi,v 2 ) + B(t; v 2 ,h 200 ) + 2B(t;vi, h lw ) - 2a 20 o A{r)v 2 )dr, 



(3.37) 

taking the normalization from (3.25) into account. Now, h 2 \Q is defined by (3.36) up 
to a multiple of v 2 , therefore we impose the orthogonality condition with the adjoint 
eigenfunction v\ to obtain h 2 i$ as the unique solution of 

/1210 - A{T)h 2W + iujh 210 - C(t; Vi, v\, v 2 ) - B(t; v 2 ,h 20 o) 

-2B{r; Vx,h no ) + 2a 200 A(T)v 2 + 2b 210 v 2 = 0, r G [0, T], 

h 2 w(T) - h 2l0 (0) = 0, 

Jo (^2^210) dr = 0. 

(3.38) 

Since is not a term in the normal form (2.4), we will find a non-singular equation 
for /1120 when collecting the ^i^|-terms from the homological equation. Moreover, because 
of property (3.22) h\ 2 Q is anti-periodic and thus 

A120 - A(r)hi 20 + 2iujhi 20 - C(r; v 1 ,v 2 ,v 2 ) 

-B(T;v 1 ,h 020 )- 2B{r;v 2 ,h ll o) = 0, re [0,7], (3.39) 
^l2o(T) + h 120 (0) = 0. 

The two remaining third order terms are the £2 (£2 1 2 -terms and the £1 |^2| 2 -terms, which 
both give a singular equation, namely 

ho2i - A(T)h<m + iwhtm = C(r; v 2 , v 2 , v 2 ) + B(t; v 2 , h 020 ) + 2B(t; v 2 , h u) 

— 2aonv 2 — 26021^2 — 2iLoaonv 2 

and 

h lu - A(r)hm = C(r; vi,v 2 ,v 2 ) + B(t; vx,h on ) + B(t; v 2 , h wl ) + B(r; v 2 , h 110 ) 
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The first function is T-periodic, the second one is anti-periodic. Both involve a Fredholm 
solvability condition, which leads to the computation of the two remaining unknown third 
order normal form coefficients of (2.4), i.e. 



1 f 

bo2i = o / ( v 2> c (t; v 2 ,v 2 ,v 2 ) + B(t; v 2 , h 020 ) + 2B(t; v 2 , h u) - 2a 01 iA(T)v 2 ) dr 
1 Jo 



(3.40) 

and 



T 

am = I (v*,C(t; vi,v 2 ,v 2 ) + B(t; vi,h 01l ) + 2^(B(t; v 2 , h 10 i)) - a Q uA(r)vi) dr. 



o 



(3.43) 



(3.41) 

Since we need both /1021 and hm for the computation of higher order normal form coeffi- 
cients, we also write down their BVPs 

h 021 - A(T)h 021 + iujh 02 i - C(r; v 2 ,v 2 ,v 2 ) - B(t; v 2 , h 020 ) 

2B(T;v 2 ,h 011 ) + 2a 0ll A{T)v 2 + 2b 021 v 2 = 0, re[0,T], 
fc<m(T)-fco2i(0) = 0, 
J (v 2 ,h 02 i) dr = 

(3.42) 

and 

A m - A(r)hm - C(r; vi, v 2 , v 2 ) - B(t; vi,h u) 

-2»(B(r; v 2 , hm)) + aonA(r)vi + amv x = 0, r£ [0, T], 

/iiii(T) + /i m (0) = 0, 
fo(vt,h in ) dr = 0. 

The stability of a possibly existing torus depends on the fourth and fifth order coeffi- 
cients, which we have listed in Appendix C. 

3.3. NSNS. The five-dimensional critical center manifold W c (r) at the NSNS bifurcation 
can be parametrized locally by (£,r) G C 2 x [0, T] as 

u = u {t) + 6^i(t) + £i«i(t) + &v 2 (t) + 1 2 v 2 {t) + H(£, r), (3.44) 

where H satisfies H(£,T) = H(^,0) and has the Taylor expansion 

H (M = E -^h m {r)^{^ 2 + Om% (3-45) 

2<i+j+k+l<5 

where the complex eigenfunctions v\ and v 2 are given by 

v\ — A{t)v\ + ioj\V\ = 0, t e [o, r], 

vi(T)- Vl (0) = 0, (3.46) 

/ {vx,v\)dr -1 = 0, 
and 

v 2 - A(t)v 2 + ioj 2 v 2 = 0, r G [0, T], 

v 2 (T)-v 2 (0) = 0, (3.47) 

J (v2,v 2 )dr -1 = 0. 

The functions v\ and ^2 exist because of Lemma 2 of [18]. The functions hij^i will 
be found by solving appropriate BVPs, assuming that (2.1) restricted to VF c (r) has the 
normal form (2.5). 

The coefficients of the normal form arise from the solvability conditions for the BVPs 
as integrals of scalar products over the interval [0, T]. Specifically, those scalar products 
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involve among other things the quadratic and cubic terms of (2.2) near the periodic so- 
lution no, the eigenfunctions v\ and v 2 , and the adjoint eigenfunctions if*, v% and v\ as 
solution of the problems 

<p*+A T (r)<p* = 0, t G [0, T], 

p*(T) -¥>*(()) = 0, (3.48) 
£ (<p*,F( Uo ))dT-l = 0, 

uj +A t (t)vI +iunvl = 0, re[0,T], 

t4(T)-^(0) = 0, (3.49) 



and 



v* 2 + A t {t)v% + iu 2 v* 2 = 0, r G [0, T], 

«2(r)-«5(0) = 0, (3.50) 

£ { V * 2 ,V 2 )dT-l = 0. 

By collecting the constant and linear terms we get the identities 

iiQ = F(u ), vi + icoivi = A(t)vx, v 2 + w 2 v 2 = A(r)v2, (3.51) 

and the complex conjugates of the last two equations. (3.51) merely reflects the definition 
of uq and the first equations in (3.46), (3.47). 

By collecting the £f (or £|-terms)-terms we find an equation for /12000 

/12000 - A(T)h 2 QQQ + 2iuih 20 QQ = B(t; vi,vi), 

(or its complex conjugate). This equation has a unique solution /12000 satisfying h 2 QQo(T) = 
^2000(0), since due to the spectral assumptions e 2tuJlT is not a multiplier of the critical cycle. 
Thus, /12000 can be found by solving 

^2000 - ^(t)^2000 + 2io;i/i 2 ooo - B(t;vi,vi) = 0, r e [0, T], , . 

^2000(^-^2000(0) = 0. K6 * Z} 

The function /10200 is just the complex conjugate of the function /i 2 ooo- Analogously, by 
comparing the £2 -terms, we find that /10020 is the unique solution of 

ft-0020 - A(T)h 0020 + 2iw 2 /ioo20 - B(t; v 2 ,v 2 ) = 0, r G [0, T], , ^ 

^0020(^-^0020(0) = 0. 

By collecting the |£i| 2 -terms we obtain a singular equation, as expected since this term 
is present in the normal form (2.5), namely 

ft-noo - ^4(t)^iioo = B(t; vi,vi) - anoo^o, 

to be solved in the space of functions satisfying hnoo(T) = /inoo(O)- Since the null-space 
is spanned by uq, the Fredholm solvability condition 

r-T 

(ip*,B(T;vi,vi) - anoouo) dr = 







gives us the possibility to calculate parameter anoo due to the normalization condition in 
(3.48), i.e. 

"1100= / (<p*,B(t;vi,vi)) dr. (3.54) 



Function /inoo is now determined up to the addition of a multiple of uq. As always, we will 
add an orthogonality condition, in this case with the adjoint eigenfunction coresponding 
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with multiplier 1. Therefore, with the value of ctnoo from (3.54) we obtain /inoo as the 
unique solution of the BVP 



Anoo - A(r)h 110 o - B(t-, vi,v{) + anoo^o 
hnoo(T) - /inoo(O) 
fo(<p*,hnao) dr 

Analogously, function /loon can be obtained by solving 

froon - A(T)h 00 n - B(t; v 2 , v 2 ) + a onu 
hoou(T) - ftoou(O) 



0011, 



dr 



0, t € [0, T] , 

0, 

0. 



0, t G [0, T], 

0, 

0. 



(3.55) 



with 



(3.56) 



(3.57) 



n lin i | - | ((p*,B(t;v2,v 2 )) dr. 
By collecting the £i£ 2 -terHis we find the following differential equation for /iioio 
/11010 - A(T)h ww + i^i/iioio + ^2^1010 = B(t; v 1 ,v 2 ). 

This equation has a unique solution /11010 satisfying /iioioCO = ^-ioio(O), since due to the 
spectral assumptions e t ( UJ i+ UJ 2) T [ s no t a multiplier of the critical cycle. Thus, /iioio can be 
found by solving 



h 



1010 



^(t)/iioio + iuJihioio + itt) 2 /iioio - S(r;ui,« 2 ) 

^ioio(r) - /iioio(O) 



0, r € [0,T], 
0. 



(3.58) 



We note that /10101 = ^1010 • 

The last second order derivative coming from looking at the £i£ 2 -terms results in a 
non-singular differential equation, such that 



hiooi - A(T)hiooi + iwi^iooi - iw 2 ^iooi - B(t;v 1 ,v 2 ) = 0, t e [0,T] 

^looi(T)-/tiooi(0) = 0. 



(3.59) 



We now come to the third order terms. From the £f and ^-terms we immediately get 
the BVPs for /13000 an d /ioo30> namely 



and 



^13000 - -4(t)^3000 + 3iwi/i 30 oo - C(t; vi,vi,vi) 

-3B(r;«i,7»200o) = 0, r£ [0,T], 
/l30O0(T)-/l3000(0) = 



^0030 - -4 (r) /i 00 30 + 3iw 2 /ioo30 - C(t;v 2 ,v 2 ,v 2 ) 

-3B(T;v 2 ,h 0020 ) = 0, r G [0,T], 
^0030 (r) - /lQ030(0) = 0. 



(3.60) 



(3.61) 



Since the £1 |£i| -term is present in the normal form for the double Neimark-S acker bifur- 
cation, a Fredholm solvability condition is involved coming from the following differential 
equation for /12100 



/12100 - A(r)h 2 ioo + w^i/12100 



C(t; vi,vi,vi) +2B(T;v 1 ,hi W0 ) + B(t; v x , /t 20 oo) 
-202100^1 - 2iwianooUi - 201100*1- 
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The differential operator -£p — A(r) + ioj\ is singular with its null-space spanned by v\, so 
we get condition 

T 

(u*,C(r; + 2B(T]Vi,h 1100 ) + B(t; «i, /i 2 ooo) - 2a 2W oVi 

-2auoo«i - 2^101100^1) dr = 0. 

Taking the normalization condition from (3.49) and the differential equation from (3.46) 
into account, we get 



1 



T 



fl2ioo = 2 / ( v *> c ( T i v i> v i> v i) + 2S(r;vi, /inoo) + B(r;«i,/i 2 ooo) - 2anoo^.(r)ui) cZt. 







(3.62) 



Therefore, we can compute /12100 as the unique solution of the BVP 



^2100 - -4(^)^2100 + ^1/12100 - C(r; ui,ui,ui) 
-2B(t; vt, /inoo) - ^(t; ui, ^2000 ) + 2a 2 ioo«i + 2aiioo^(r)ui = 0, r G [0, T], 

folQoCD - /l2100(0) = 0, 

/ T K,/i 2 ioo) dr = 0. 

(3.63) 

We can now immediately list the following four BVPs 

/12010 - -4(r)/i 2 oio + 2iwi/i 2 oio + ^2^2010 - C(r; v 1 ,vi,v 2 ) 

-B{r;v 2 ,h 2000 )- 2B{T;v l ,h l0l0 ) = 0, rG[0,T], (3.64) 
^2010 (r) - ^2010(0) = 0, 

^2001 - ^(t)^2001 + 2io;i/i 2 ooi - iw 2 /i 2 ooi - C(r; u 1; «i, « 2 ) 

-B(r;w2,/i20oo)-2B(r;wi,/i 10 oi) = 0, r e [0,T], (3.65) 

^2001 (r) - /l200l(0) = 0, 

hiow - A(r)fi,io2o + iwihiom + 2iuj 2 h 1020 
-C{r;v l ,v 2 ,v 2 )-B(T;v 1 ,h 0020 )- 2B(T;v 2 ,h 1010 ) = 0, t € [0,T\, (3.66) 

^102o(r) - /il02o(0) = 0, 
and 

^0120 - -4(t)^oi20 - ^1^0120 + 2iu 2 h i 20 
-C(t;vi,v 2 ,v 2 ) - B(r]V 1 ,hoa2o) - 2B(T;v 2 ,h 0110 ) = 0, rG[0,T], (3.67) 

/i i2o(r) - ^0120(0) = 0. 

The £ 2 |£ 2 | 2 -terms from the homological equation make it possible to compute 60021- 
Indeed, the differential equation 

^0021 - A(r)hoo2i + iw 2 /ioo2i = C(r; v 2 , v 2 , v 2 ) + B(t; v 2 , h 0020 ) + 2B(t; v 2 , h 00 n) 

-2b 002 iv 2 - 2aooii^2 - 2iw 2 aooii^2 

results in a solvability condition with v 2 , i.e. 

/ (v 2 ,C(t;v 2 ,v 2 ,v 2 ) + B(T;v 2 ,h 0020 ) + 2B(t; v 2 , h on) 
Jo 

-2b o 2 \v 2 - 200011^2 - 2iuj 2 a onv 2 ) dr = 0. 
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Therefore, considering the normalization condition from (3.50) and the differential equa- 
tion from (3.47), we can calculate parameter 60021 



1 



T 



60021 = - / (v2,C(t;v2,v 2 ,v 2 ) + B(t;v 2 , /10020) + 2-B(r; v 2 , h ou) - 2a 00 iiA(T)v 2 ) dr, 







(3.68) 

with /10021 as the unique solution of the BVP 

%)2i - A(T)h 0021 + iuj 2 h 0021 - C(t; v 2 ,v 2 ,v 2 ) 
-B(T;v 2 ,h 0020 ) - 2B(T;v 2 ,h 0011 ) 

+2b 002 iv 2 + 2aooiiA(T)v 2 = 0, r e [0,T], (3.69) 
hoo2i(T) - /i 02i(0) = 0, 
Jo(v 2 ,h 002 i)dT = 0. 

The last two third order terms which we have to examine give us both the formula for a 
normal form coefficient. The first one, obtained from the |£i| 2 ^2-terms, gives us the BVP 

/lino - A(T)h 1U0 + iuj 2 h 1110 - C(r; vi,vi,v 2 ) - B(t; vt, fyjlio) 

-B(t; v 1} h ww ) - B(t; v 2 , /inoo) 

+61110^2 + ano A(t)v 2 = 0, t£[0,T], 
h n w(T) -/imo(0) = 0, 
Jo (v 2 , h mo) dr = 0, 

(3.70) 

where from the solvability condition it follows that 



mio = / (v 2 ,C(t;vi,vi,v 2 ) + B(t; v 1} /toiio) + B(t; vi, frioio) + B(t; v 2 , /inoo) - o>i 10 oA(t)v 2 ) dr. 
Jo 



Analogously, we obtain the BVP 

frion - A(r)h wll + iuJthwu ~ C(r; v 1 ,v 2 ,v 2 ) - B(t; vi,h 00 n) 

-B(t; v 2 , /iiooi) - B(r; v 2 , h 10W ) + aion^i + a oii-4(T")ui = 0, r G [0, T], 

hiou(T) - /iion(0) = 0, 

Jo ( v h h wn) dr = 0, 

(3.71) 

with 



■T 

aion = / (vI,C(t; vi,v 2 ,v 2 ) + B(t; v x , h 00 u) + B(t; v 2 , /iiooi) + B(r; v 2 , /iioio) - 0£ oiiA(r)vi) dr. 



o 



We still need the coefficients 6noi and aom which are determined by 

fT 



?1101 



/ (v 2 ,C(t;vi,vi,v 2 ) + B(t; vx,h ioi) + B(T;vi,h 100 i) + -B(r; v 2 , /inoo) 
Jo 

-a noo A(r)v 2 ) dr (3.72) 



and 



aom = / (vI,C(t;vi,v 2 ,v 2 ) + B(t; «i, /loon) + B(r;v 2 ,h ioi) + B(t; v 2 , /i no) 
Jo 

-Q!ooii^(t)^i) dr. (3.73) 

As before, the higher order terms which determine the stability of the torus can be 
found in Appendix C. 
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3.4. Implementation. Numerical implementation of the formulas derived in the previous 
section requires the evaluation of integrals of scalar functions over [0, T] and the solution 
of nonsingular linear BVPs with integral constraints. Such tasks can be carried out within 
the standard continuation software such as AUTO [11], content [25], and matcont [8]. 
In these software packages, periodic solutions to (1.1) are computed with the method of 
orthogonal collocation with piecewise polynomials applied to properly formulated BVPs 
[0,4]. 

We have implemented our algorithms in matcont analogously to the eight cases with 
n c < 3. For further details we refer to [7] where this is extensively discussed. 

4. Examples 

4.1. Laser model. In [34] a single- mode inversionless laser with a three- level phaser was 
studied and shown to operate in various modes. These modes are "off' (non-lasing), con- 
tinuous waves, periodic, quasi-periodic and chaotic lasing. The model is a 9-dimensional 
system given by 3 real and 3 complex equations: 

Paa =Ra~ |(^K& ~ Kb) + %{°ac ~ Kc)) 
< Pbb = Rb + ^l(0- a b ~ 0-* ab ) ^ 
O-ab = -(71 + iAl)(Tab ~ %(fll(paa ~ Pbb) ~ QpVcb) 
<?ac = -(72 + iA p )fJ ac - ^{Vt p {2p aa + p bb - 1) - 
Ocb = -(73 + ~ A p ))<T cb - §(0/(7* c - %0~ab), 

with R a = -0.505p aa - 0A05p bb + 0.45, i? 6 = 0.0495/> aa - 0.0505pm, + 0.0055 and A* = 
&-cav + gyi{o~ab)Qh The fixed parameters are 71 = 0.275,72 = 0.25525,73 = 0.25025, ^ cav = 
0.03,5 = 100, A p = 0. The parameters Q p and A cav are varied. The bifurcation diagram 
of (4.1) is computed in [28] and is reproduced in Figure 5. 

4.1.1. The LPNS points. Figure 5 shows three NS curves NS(1), NS(2) and NS(3) start- 
ing from two HH points. On NS(3) one of the richer situations happens. The normal 
form coefficients for the LPNS point at (£l p ,A cav ) = (3.411,-1.819) are (s,0,E) = 
(1,-0.139,-911.248), so s9 < 0. This means that there exists a 3-torus, which is sta- 
ble since 9 < and E < 0. Therefore, we are in the case represented in Figure 1 (c), but 
with a stable 3-torus. For computing the Lyapunov exponents, we used a code written by 
V. N. Govorukhin (2004). Figure 6 (left) shows the calculated Lyapunov exponents for 
Q p fixed at 3.45 and A cav £ [—1.8; —1.6]. More detail is shown in Figure 6 (right), where 
we get a clear view on the number of Lyapunov exponents equal to zero. For A cav values 
to the right of —1.636, there is one Lyapunov exponent equal to zero, which corresponds 
to the stable limit cycle from region 6 in Figure 1 (c). At A cav = —1.636, we cross NS(3) 
and arrive in region 5 with a stable 2-torus and therefore two Lyapunov exponents equal 
to zero. When crossing the P curve at A cav = —1.773, the stable 3-torus from region 4 
arises. Remark that in some small intervals only two Lyapunov exponents are equal to 
zero, and thus not the expected three zero ones, but these correspond with resonances on 
the 3-torus. Then, in the interval A cav G [—1.796; —1.7916] positive Lyapunov exponents 
appear which indicates that there is chaos. This zone corresponds with T. Afterwards, 
we arrive in region 3, where all Lyapunov exponents are smaller than zero. 

On the NS(2) curve there is one LPNS point for (Cl p , A cav ) = (4.632, 1.438). The normal 
form coefficients are (s,6,E) = (1,0.206,808.009). The product s9 > is positive, so we 
are in a "simple" case, where no 3-torus is present. Since s = 1, the torus arisen through 
the Neimark-S acker curve exists below the NS(2) curve. We have computed the Lyapunov 
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Figure 5. Bifurcation diagram of (4.1). The thin red curves are Hopf 
curves. In blue are limit point of cycles bifurcations and in magenta 
Neimark-Sacker bifurcations. Solid/dotted curves correspond to supercrit- 
ical/subcritical bifurcations. The dashed curves are curves of neutral sad- 
dles. 



0.01 




1.785 



Figure 6. Lyapunov exponents computed for Q p = 3.45 close to the LPNS 
point at (n p ,A cav ) = (3.411,-1.819), (left) for A cav G [-1.8; -1.6] and 
(right) zoomed in near the region with chaos due to heteroclinic tangles. 
The vertical black lines indicate the parameter values where a bifurcation 
occurs. 
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exponents for a straight line where the beginning point (tt p ,A cav ) = (4.302,0.673) and 
end point (Q p ,A cav ) = (4.984,1.984) lie between the curves LPC(2) and NS(2), to the 
left and to the right of the LPNS point. In Figure 7, we plot the Lyapunov exponents 
for ftp E [4.3,4.98]. The stable limit cycle is situated in the upper wedge between the 
LPC(2) and NS(2) curves which corresponds to region 4 in Figure 1(a), so we have one 
Lyapunov exponent equal to zero for Q p -values larger than the subcritical NS(2) curve. 
At £l p ~ 4.41, we cross the subcritical NS(2) curve, with to the left no zero Lyapunov 
exponents. 



0.02 






4.5 



4.7 



Q 



Figure 7. Lyapunov exponents computed close to the LPNS point at 
(Q p ,A cav ) = (4.632,1.438). The two-coloured dashed lines reveal pairs of 
equally large Lyapunov exponents. 



4.2. Periodic predator-prey model. As a second model we study a simple two-patch 
predator-prey system with periodic (seasonal) forcing. Simple predator-prey models lead 
to the 'paradox of enrichment', i.e., increasing the carrying capacity of the prey ultimately 
leads to extinction of the population [30]. Outside the laboratory, however, stable pop- 
ulations are observed and not extinction. Here, spatial models have been put forward 
to explain this discrepancy. As the simplest spatial case, one may consider a two-patch 
predator-prey model [20] where predator and prey can migrate between the two patches by 
diffusion. This leads to a diffusive instability of large oscillations and stabilizes the total 
population size [21]. Here, we propose an extension where one of the patches experiences 
seasonal influences while the other can be seen as a wild-life refuge where human interven- 
tion minimizes seasonal influences. As a simplication we will only consider the case that 
the predators can move between the patches, i.e., they can cross the refuge barrier. On a 
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proper time scale, the investigated system is defined by 




2/2 = -2/2 H —r- + 1(X2 ~ V2) 

yi + b 2 



(4.2) 



Vl = — «2 + «l(l — «i — «§) 
<y 2 = Vi + V 2 (l ~v\- vl). 



The values of x\ and £2 denote the numbers of individuals (or densities) respectively of 
prey and predator populations living outside the refuge and y\ and y% are the corresponding 
numbers or densities inside. The intrinsic growth rates r{ and the constant attack rate c 
are parameters of the model. For the predator outside the refuge, the Holling type II is 
chosen as functional response with a half saturation which varies periodically with period 
2ir. To this end, the last two equations are introduced; their solutions converge to a stable 
limit cycle v\(t) = cos(t + 0) with a phase shift <p depending on the initial conditions. The 
terms with parameter 7 describe the coupling of the two patches. The fixed parameter 
values are r\ = l,r 2 = l,b% = 0.4,7 = 0.1, c = 2. We will use the half saturation b 2 as 
a continuation parameter together with the amplitude of the seasonal forcing e. It is not 
our aim to give a full analysis of this model, but rather analyze the codim 2 bifurcations 
relevant for this paper. We observe that a refuge can induce complex behaviour in a spatial 
population model with seasonal forcing. 

4.2.1. The P DNS points. Figure 8 represents a bifurcation diagram for system (4.2) where 
two PDNS points are detected. The right PDNS point has parameter values (62,6) = 
(0.277, 0.530). We are in the "simple" case of Section 2.2.2 because the product of the co- 
efficients pu = —5.01 ■ 10 -2 and P22 = —0.211 is positive. Since 9 = —0.320 and 5 = 1.087, 
Figure 2(a) indicates that the bifurcation diagram in a neighbourhood of the PDNS point 
is as in case III in Figure 3(a), where [i\ = corresponds with NS1 and H2 = with PD. 
Curve T\ corresponds to the Neimark-S acker curve of the period doubled cycle NS2(2) in 
Figure 8. Therefore, we expect the period doubling 'curve' T2 of the torus to be situated 
to the left of NS1(2) and under the PD curve. The stable limit cycles are situated in the 
lower right region of the PDNS point. The exact location of T2 can be determined by 
computing Lyapunov exponents for fixed b 2 values smaller than the critical 62 = 0.277 
corresponding with the PDNS point. We have plotted a sketch of this T2 curve in Figure 

9, which represents a zoom of the neighbourhood of the PDNS point and which includes 
a plot of NS2(2) (curve T\ in Figure 3(a)). We have computed the Lyapunov exponents 
for 62 fixed at 0.261 and e £ [0.46; 0.62], see Figure 10. In this figure the black vertical 
lines indicate the position of the PD and NS2(2) curves. From the value of the Lyapunov 
exponents we derive that T2 is crossed for e ~ 0.52. To the left of the T2 curve in Figure 

10, we have a stable torus, arisen through the supercritical Neimark-S acker curve NS1(2), 
corresponding with region 2 from Figure 3(b). Between the curves T2 and NS2(2), the 
2-torus arisen through T2 is attracting. These regions correspond with region 6 (between 
T2 and PD) and 5 (between PD and NS2(2)) from Figure 3(b). When crossing the NS2(2) 
curve, the 2-torus disappears and the period doubled cycle becomes attracting. All this is 
in agreement with the fact that two Lyapunov exponents are equal to zero to the left of 
NS2(2), where afterwards only one zero Lyapunov exponent is left. 



26 



V. DE WITTE, W. GOVAERTS, YU.A. KUZNETSOV, AND H. G. E. MEIJER 




0.2 - 
0.1 - 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 

b 2 

Figure 8. Bifurcation diagram of limit cycles in (4.2). In green are period 
doubling curves and in magenta Neimark-S acker curves (of the first or of 
the second iterate, respectively labeled with NS1 and NS2). 

The left PDNS point at (b 2 ,e) = (8.699 • 1CT 2 , 0.519) again belongs to one of the "sim- 
ple" situations in Section 2.2.2 (pn = — 0.447, P22 = —1.472). The neighbourhood of the 
bifurcation point is as in case I in Figure 2(a) since (6,5) = (2.234,1.304). Remark that 
the stable limit cycles are situated in the lower left quadrant of the PDNS point in Figure 
11. The behaviour in a neighbourhood of this PDNS point can be derived from Figure 
11, which includes a plot of the Neimark-S acker curve NS2(1) of the period doubled cycle 
and also a sketch of the period doubled curve T2 of the torus, made on the basis of the 
computation of the Lyapunov exponents. We have calculated the Lyapunov exponents for 
parameter values in the upper right quadrant, close to the PDNS point, for 62 = 0.08709. 
The results are given in Figure 12. Going from the left to the right, where we follow the 
solid lines, we start with two Lyapunov exponents equal to zero which correspond with 
the stable torus from the original cycle in the regions 2, 3 and 4 from Figure 3. At the 
point where the second Lyapunov exponent becomes non-zero, the T2 curve is located, 
namely at e ~ 0.5198. We then arrive in region 12 from Figure 3(b) where the 2-torus 
has lost his stability and the period doubled cycle is stable. Therefore, one zero Lyapunov 
exponent remains. We scan the Lyapunov exponents for a second time where we now go 
from the right to the left and follow the dashed lines. The second Lyapunov exponent now 
approaches zero not at the T2 curve but at the NS2(1) curve. This is explained by the 
bistability happening in region 4, where one Lyapunov exponent equal to zero indicates the 
stable period doubled cycle and two zero Lyapunov exponents indicate the stable torus. 
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Figure 9. Zoom of the neighbourhood of the PDNS point at (62,6) = 
(0.277,0.530) from Figure 8. In green are period doubling curves, in ma- 
genta Neimark-S acker curves (of the first or of the second iterate, respec- 
tively labeled with NS1(2) and NS2(2)), in blue is the sketch of the T2 
'curve'. 
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Figure 10. Lyapunov exponents computed for b 2 = 0.261, close to the 
PDNS point at (b 2 ,e) = (0.277,0.530). 



2N 



V. DE WITTE, W. GOVAERTS, YU.A. KUZNETSOV, AND H. G. E. MEIJER 




Figure 11. Zoom of the neighbourhood of the PDNS point at (62, e) = 
(8.699 • 10 -2 , 0.519) from Figure 8. In green are period doubling curves, 
in magenta Neimark-S acker curves (of the first or of the second iterate, 
respectively labeled with NS1(1) and NS2(1)), in blue is the sketch of the 
T2 'curve'. 

When going further, we cross region 3 and 2, with the stable torus of the orginal cycle. 

Remark that since we have a periodically forced system the return time is independent 
of the distance from the limit cycle, so we could do this extra check. Indeed, for all 
PDNS points, the anjk in the first equation of (2.4) are zero up to the accuracy of the 
computation. Here too, the Lyapunov exponents corroborate the prediction based on the 
normal form coefficients. 

4.3. Control of vibrations. In [ ] a two-mass system of which the main mass is excited 
by a flow-induced, self excited force is studied. A single mass which acts as a dynamic 
absorber is attached to the main mass and, by varying the stiffness between the main mass 
and the absorber mass, represents a parametric excitation. The system is given by 

xi = Vl 
£2 = v 2 

vi = -h(vi - v 2 ) - Q 2 (l + eyi)(xi - x 2 ) , 
* v 2 = Mki(vi - v 2 ) + MQ 2 (1 + eyi)(xi - x 2 ) - k 2 v 2 - x 2 + f3V 2 (l - jv%)v 2 
Vl = ~VV2 + - y\ - yl) 
= mix +2/2(1 -yj- yl). 

The following parameters are fixed: e = 0.1, = 0.1, /3 = 0.1, V = V2A, 7 = 4, Q = 
0.95, M = 0.2, ki and ij will be the continuation parameters. 
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Figure 12. Lyapunov exponents computed for b 2 = 0.08709, close to the 
PDNS point at (b 2 ,e) = (8.699-10" 2 , 0.519). Exponents indicated with solid 
lines are computed by following the attractor with increasing e, dotted lines 
with decreasing e. This highlights the bistability between NS2(1) and T2. 

4.3.1. The NSNS points. An NSNS point is detected for (ki,r)) = (9.167 • 10" 2 , 0.411), see 
Figure 13. The normal form coefficients are 




Figure 13. Partial bifurcation diagram of limit cycles in system (4.3). 



(pu,P22, 0, 5, sign h) = (-3.733 • 10~ 3 , -6.494 • 10" 3 , 0.541, 1.203, 1). 

The positive sign of the product P11P22 implies that we are in a "simple" case of Section 
2.2.3. Since 5 > 9, the role of both coefficients has to be reversed. Therefore, 6 > 1,6 < 
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1,95 < 1 indicate that the NSNS bifurcation is located in region II in Figure 3(a). As in 
the previous examples, we have computed the Lyapunov exponents to check the obtained 
results of the normal form coefficients. We have done the computations for k\ fixed at 0.083 
and T) G [0.4; 0.42] [r\ values are between the NS curves). The results are given in Figure 
14. For jy-values starting from 0.38, we are in region 3 (or 12 due to symmetry) in Figure 
3(b), where there is a stable 2-torus and thus two Lyapunov exponents equal to zero. A 
third Lyapunov exponent approaches zero and between ry ~ 0.4117 and ry ~ 0.4154 three 
Lyapunov exponents are equal to zero. This region denotes the appearance of a stable 3- 
torus and corresponds with region 5 from Figure 3(b) II. The critical values of rj correspond 
with the curves T\ and T<i in Figure 3(a). For rj > 0.4154, only a stable 2-torus remains 
such that there are two zero Lyapunov exponents. Therefore, the computed Lyapunov 
exponents are in agreement with the normal form coefficients. 




0.4 0.405 0.41 0.415 0.42 



Figure 14. Lyapunov exponents computed for k\ = 0.083. 

Also in this case all oiijki in the normal form (2.5) vanish since we have a periodically 
forced system. 

5. Discussion 

This paper completes the development of efficient methods for the computation of the 
critical normal form coefficients for all codim 1 and 2 local bifurcations of limit cycles, 
started in [24, 7] and based on [18]. Together with our previous papers on the computation 
of the critical normal form coeffcients for codim 1 and 2 local bifurcations of equilibria in 
ODEs [22] and fixed points of maps [26, 27], it contributes to the development of methods, 
algorithms, and software tools for multiparameter bifurcation analysis of smooth finite- 
dimensional dynamical systems. 

The resulting formulas are independent of the phase space dimension and can be ap- 
plied in the original basis, without preliminary linear transformations. As limit cycles 
are concerned, the formulas are directly suitable for numerical implementation using or- 
thogonal collocation. They fit perfectly into a continuation context, where limit cycles 
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and their bifurcations are computed using the BVP-approach [12], without numerical ap- 
proximation of the Poincare map or its derivatives. Being implemented into the matlab 
toolbox matcont [8, 9], the developed methods are freely available to assist an advanced 
two-parameter bifurcation analysis of dynamical systems generated by ODEs and maps 
from various applications. 

To fully support the two-parameter bifurcation analysis of ODEs and maps, one needs 
special methods to switch between various branches of codim 1 bifurcations of fixed points 
and cycles rooted at codim 2 points. Such methods have been developed and implemented 
in matcont for codim 2 equilibrium [28] and fixed point [14] bifurcations. Switching at 
codim 2 points to the continuation of codim 1 local bifurcations of limit cycles seems to 
be the next natural problem to attack, while that for codim 1 bifurcations of homoclinic 
and heteroclinic orbits is more difficult and probably requires new ideas. Similar remarks 
can be made about quasiperiodic bifurcations of tori. 



Appendix A. Derivation of the critical normal forms 

A.l. Notation. Let M G M nxri be the monodromy matrix. In all codimension 2 cases 
all critical multipliers, i.e., all multipliers with modulus 1, have non-degenerate Jordan 
blocks. Let Mq be the critical Jordan structure, i.e., the block diagonal matrix consisting 
of the critical Jordan blocks, starting with the block of the trivial multiplier 1. Let [1% = 
e iek (0 < 9k < 2tt) be a critical multiplier with multiplicity rrik- The matrix E 



is defined as 







1 



0\ 




: •. 1 

\0 ... a k ) 

where is the Floquet exponent of the multiplier fj,^, with = iO^/T in the case of a 
positive real multiplier or a complex multiplier and = for fi^ = — 1. The matrix Lo 
is the block diagonal matrix formed from the blocks for which \fj,^\ = 1, starting with 
the block that corresponds with multiplier 1. The matrix Lq is the matrix Lq without the 
first row and the first column. 



A. 1.1. LPNS. At the LPNS bifurcation the matrices described in A.l are 



M n 



A 1 

1 

e 

\0 






iuoT 






-iuT J 



u 





1 





^ 




















iuj 













— ioj J 







-iu i 



We are in a case in which we can apply Theorem 2 from [18]. So we can define a T-periodic 
normal form 

where £ = (£i, £2, ^2) T with £1 6l,^ £ C. The polynomials p and P are real, respectively 
complex, T-periodic in r and at least quadratic in (£i,£2>£2) such that 

^P(r, - ^P(r, C)L^ = 0, |-P(r, fl + L* P(r, - -|p(r, fl^S* = 0. 
If we write the polynomials in a Fourier expansion, namely 

oo oo 



2-kIt 

e" t . 
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we obtain for any I £ Z the following differential equations 

± m) iii-t^m)-Lim) = v. 

Putting Lq into the equations and writing 

P(6,6,6) = (^ 1) (ei,6,6), J P / (2) (6,6,6),p/ 2) (ei,6,6)) T 

we can rewrite them as a set of differential equations in variable £2 

* w 6-^pK6,6,6) + «¥m6,6,6) = ^6-J^(6,6,6), 

d6 T d& 

ico& A P (i) (^,6) + »^ (1) (6, 6,6) = ^ 2 -J^ (1) (6,6,6), 
«6 J «6 

iw6^-P, (2) (6,6,6) + i^ (2) (6,6,6) = ^6-^ (2) (6,6,6) + ^ (2) (6,6,6), 

^6 ^-p/ 2) (6,6,6) + ^p/ 2) (6, 6, 6) = ^6-iif >(6, 6, 6) - ^p/ 2) (6, 6, 6)- 
«6 J a6 

k(6,6,6),p/^(6,6,6) and P/^(6,6,6) are polynomials, and from the equations it 

follows that w(6,6,6),P/ 1) (6,6,6) and P/ 2) (6,6,6) are zero if I ^ 0. Thus the 
polynomials are r-independent. We obtain 

6-j7-po(6,6,6) = 6^po(6,6,6), 
«6 "6 

6^7-- p o^(6,6,6) = 63^-^0^(6,6,6), 
«6 "4 2 

6j- p o (2) (6,6,6) = 6^-p (2) (6,6,6) + p (2) (6,6,6), 

"6 «6 
and the complex conjugate of the last equation. From the first equation it follows that 

po(6,6,6) = Y>i(6) + MI6I 2 ) + ^3(6)MI6I 2 ), 

where ip2, "03 and -04 are at least linear in their argument and ipi at least quadratic. 
Similarly, we obtain 

Po (1) (6,6,6) = 0i(6) + ^(|6I 2 ) + ^3(6)04(|6I 2 ), 

with the same conditions for <j) as the ones for ip. At last, from the third equation we can 
derive that 

P (2) (6,6,6) = 6x1(6) + 6x 2 (|6l 2 ) + 6x 3 (6)x4(|6l 2 ), 

where Xi,X2,X3 and X4 are at least linear in their argument. 

Assembling all the information gives us the following normal form 

— = 1 + 6 + «2oo6 2 + "on I6I 2 + 03006* + «'ni6 16I 2 + • • • , 

= a 200?l + a 'oil I6I 2 + O3006 3 + Olll6 I6I 2 + • • • , 

^p- = iw6 + &no66 + &210C16 + &0216 16I 2 + • • • • 
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We do the substitution £1 \-t — £1 and find the Iooss normal form (2.3), i.e. 

^ = 1 - 6 + a 2 oo£i + "oil I6| 2 + 0300^1 + "1116 |6| 2 + • • • , 
= ^200^1 + a on |6| 2 + ^300^1 + am6 |6| 2 + • • • , 
^ ^7 = + &11066 + &2io£i6 + &0216 16 1 2 + ■ ■ ■ , 

where the dots denote 0(|£| 4 ) terms. Note that the time evolution can be obtained by 
applying the chain rule to this system and is then given by 

' dr_ 
~dt 



1 -£l + "200£l +"011 |6| +«30oCl + «Hl6 161 +•••, 



-ijj- = 0200^1 + «on I6I 2 + a' 3 ooCi + a'in6 |6| 2 + • • • , 



with a'n. 



dt 

-0200 + a 30 o and a' m 



^ ^ = + &11066 + &2io£l6 + &0216 I6I 2 + ■ ■ ■ 

sui, and with b\ 



110 



-iuj 



&110,&210 



^300 — — "200 T "300 "HI — — O011 

iwa^oo — &110 + ^210 and b' 021 = iwaon + &021- We could use this system as our starting 
normal form. To draw conclusions about the bifurcation diagrams, we could then perform 
the time reparametrization on the center manifold to obtain an autonomous truncated 
ODE which approximates the Poincare map as a T-shift (as done in [7]). It could be 
studied by comparing the obtained ODE to the one for the Zero-Hopf bifurcation of 
equilibria. However, the time reparametrized ODE has exactly the same form as (2.3), so 
we can as well use (2.3). The approach we will follow in this paper (with (2.3) as starting 
normal form) is mathematically equivalent to the one of [7], but takes a shorter path. 



A.1.2. PDNS. At the PDNS bifurcation the matrices described in A.l are 
/1 \ /0 000\ 



Mr 







-1 
e 
\0 





iuiT 








-iu)T 










iuj 
j \o -iuj 

We are in a case in which we can apply Theorem 3 from [18]. So we can define a 2T-periodic 
normal form 

where £ = (£1, 6) 6)- The polynomials p and P are 2T-periodic in r and at least quadratic 
in their argument such that 



p{t + t, 6,6,6) 
p«(t + t, -6,6,6) 
p( 2 )( r + r,-£i,£ 2 ,e 2 ) 

and the complex conjugate of the last equation. 



0, 



p{t, -6,6,6), 

- j p (1) (t,6,6,6) 

^ (2) (r,6,6,6), 
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As in the LPNS case (since the Lq matrix is the same) we obtain that all polynomials 
are independent from r, and thus we can rewrite the last three equations as 

p(6,6,6) = p(-6,6,6), 

i> (1) (-6,6,6) = -^ (1) (6,6,6), p (2) (-6,6,6) = P (2) (6,6,6), 

thus p and P^ are even in £1 and P^ is odd in £1. Taking the results from the LPNS 
case into account, we obtain 

Potii,b,b) = Mg) + M\&\ 2 ) + Mg)M\b\ 2 ), 
^ (1) (6,6,6) = diMtf) + 6-MI6I 2 ) + 6<fe(6 2 )<MI6l 2 ), 
^o (2) (6,6,6) = 6xi(^i) + 6x 2 (|6l 2 ) + 6x 3 (6 2 )x4(|6l 2 ), 

with all functions at least linear in their argument. 

Assembling all the information gives us the Iooss normal form (2.4), i.e. 

— = 1 + 02006 + "on 161 +04006 +"022 
= a 30o£i + aiu6 I6I 2 + asooCi + 01226 |6| 4 + 0311C1 I6I 2 + • • • , 

^p- = iu& + 6210C16 + &0216 I6I 2 + &4io£ 4 6 + 6221C16 I6I 2 + &0326 16I 4 + • • ■ , 

k. or 

where the dots denote 0(|^| 6 ) terms. Note that the time evolution can be obtained by 
applying the chain rule to this system and is given by 

dr 

— = 1 + a 2 00£ 2 + "Oil 161 +«400Cl +«022 |6I +"21lC 2 |6l +•••, 

^ = asoo^i + 01116 I6I 2 + «50o6 5 + a 3 n6 3 I6| 2 + ^1226 I6| 4 + • • • , 

I H = iUJ ^ + b ' 210 ^ 2 + b ' 02lC2 16,2 + b ' 410 ^ 2 + b ' 221 ^ 2 16,2 + b ' 032 ^ 2 16,4 + ■ ■ ■ ' 

with a' 500 = 0300^200 + «500) a 3n = 0300^011 + aill"200 + «311, a'i22 = a lll a 011 + «122, 
b'210 = ^"200 + &210, b' 021 = iuaoil + &021, &410 = ^"400 + &210«200 + hio, &221 = iu}a 211 + 

&2ioaon + &02i«200 + ^221 and b' 032 = iu)a 022 + fyraaon + &032- 



A. 1.3. NSNS. At the NSNS bifurcation the matrices described in A.l are 



/I 











\ 





e iu)\T 

















e -iwiT 






















\o 











e~ iW2T ) 














\ 







































iuj 2 





v° 











-iu) 2 ) 



Lin 











^ 





— iui 














10J2 





v° 








-iu) 2 ) 



We are in a case in which we can apply Theorem 1 from [] 
normal form 

dr df 
" 1+P(r,0, -^ = LoC + P(r,6, 



So we can define a T-periodic 



dt 



d.T 
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where £ = (£i,6,6,6)- ^ ne polynomials p and P are T-periodic in r and at least 
quadratic in (£1, £1, £2, £2) such that 

£p(r, £) - |p(r, £)LS£ = 0, ^P(r, £) + L* Q P(r, £) - -|p(r, £)L^£ = 0. 

Writing down the polynomials in a Fourier expansion results in the following equations 

^16-77-^(6,6,6,6) "^26-7^(6,6,6,6) + «^rW(£l, 6, £2, £2) 

d£i a£ 2 J 

= ^i£i- r p i (£i,£i,£ 2 ,£2) + ^26-7^^(6,6,6,6), 
d£i d£ 2 

a£i d£ 2 1 

= + ^-ip/%, £1, £2, £2) + ^iP, {1) (£1, £1, £2, £2), 

d£l cf£ 2 

iwi^i ^-P/ I} (£i, £1, £ 2 ,£ 2 ) + ^ 2 £ 2 ^-p/ 1) (£i, £1, £2, £2) + iuiP^ (£1, £1, £2, £2) 

rf£l Gf£ 2 

= ^ 1 £ 1 ^P/ 1) (£ 1 ,£ 1 ,£ 2 ,£ 2 ) + za; 2 £ 2 ^P/ 1) (£ 1 ,£ 1 ,£ 2 ,^ 

d£i cf£ 2 J 

iw^i J-p/ 2) (£1, £1, £2, £2) + ^ 2 £ 2 J-p/ 2) (£i, £1, £2, £2) + ^p/ 2) (£i, £1, £2, £2) 

ef£i cf£ 2 1 

= iwi£-ip J (2) (£i,£i,6,£2) + ^ 2 £ 2 -ip/ 2) (£i,£i,£ 2 ,£ 2 ) + ^ 2 P/ 2) (£ 1 ,£ 1 ,£ 2 ,£ 2 ) 
«£i "£2 

^4^(6,(1,6,6) + ^ 2 £ 2 J-p/ 2) (£i, £1, £2, £2) + iwiP/ 2) (£1, £1, £2, £2) 
«£i "£2 

= ^^Ip^fe, £1, £2, £2) + ^ 2 £ 2 ^1p/ 2) (£i, £1, £2, £2) - ^p/ 2) (£i, £1, £2, £2). 
a£i a£ 2 i 

pK6,6,6,6),^ (1) (6,6,6,6) and p/ 2) (6,£i,6,6) are polynomials, and from the 
equations follows that / has to be equal to zero, thus the polynomials are r-independent. 
We obtain 

6-jt-po(6,6,6,6) = 6-7=^0(6,6,6,6), 
«£i a£i 

d - . - d 

6-77-M6,6,6,6) = 6-7^0(6,6,6,6), 

a£2 Ct£2 

6-J-^ (1) (6,6,6,6) = £i^Po (1) (£i,£i,£2,£2) + P (1) (£i,£i,£2,£ 2 ), 

a£i d£i 

£2-7^^0^(6,6,6,6) = 6-t^Po^(6,6,6,6), 

«?2 a£2 

6^^0^(6,6,6,6) = 6-j^Po^(6,6,6,6), 

a£i c(£i 

6-J-if ) (6,6,6,6) = £2^p (2) (£i,£i,£2,£2) + Po (2) (£i,£i,£2,£ 2 ). 

a£2 a£ 2 
From the first two equations follows that 

po(£i,£i,£2,6) = V'i(I£i| 2 ) + V'2(|£2| 2 ) + V'3(I£i| 2 )V'4(|£2| 2 ). 
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From the third and fourth equation, we obtain 

P o (1) (6,6,6,6) = 6^(l6| 2 ) + ei02(|6| 2 ) + 603(l6l 2 )</»4(|6| 2 ), 

and analogously 

Po (2) (6,6,e2,e2) = 6xi(ie 2 | 2 ) + 6x2(|6l 2 ) + 6x3(l6l 2 )x4(|6l 2 ), 

where all functions are at least linear in their argument. 

Assembling all the information gives us the Iooss normal form (2.5), i.e. 

^ -j£ = 1 + "1100 |£l| 2 + "0011 |6| 2 + "2200 |£l| 4 + "0022 |6| 4 + "1111 |6| 2 |6| 2 + • • • , 
^^sl 2 2 4 4 2 2 

— = iuiti + a 2 iooCi 161 + Oion6l6l +032006 161 +010226 161 +021116 161 161 +•••, 

— - = io; 2 6 + 600216 16I 2 + 011106 16I 2 + ^00326 16I 4 + 022106 16I 4 + 611216 16I 2 16I 2 + • • • > 

< CLT 

where the dots denote 0(|£| 6 ) terms. Note that the time evolution can be obtained by 
applying the chain rule to this system and is of the form 

2 2 ) 4 2 2 

— = 1 + Qll00|6l +"0011 161 +"2200 |6l + "0022 |6l +«llll|6l 161 + •••, 

= iw i6 + o 2100 6 16I 2 + oi u6 |6| 2 + 032006 I6I 4 + oi 2 2 6 16I 4 + o 2 m6 I6| 2 16| 2 + • • • , 
y-^r = iuj 2& + 600216 16I 2 + 6in 6 16I 2 + 600326 16I 4 + 6 2210 6 16I 4 + &mi6 16I 2 16I 2 + 

where the coefficients with primes are functions of the original coefficients. 



Appendix B. Bifurcations of the amplitude system for Hopf-Hopf 

BIFURCATION IN THE "DIFFICULT" CASE 

Here we derive quadratic approximations of the Hopf and heteroclinic bifurcation curves 
for the Hopf-Hopf truncated amplitude system (2.9), that can be written in the rescaled 
form 

x \ _ / x(/ix + x - 9y + Qy 2 ) 
y J ~ V v(t*2 + 5x -y + Ax 2 ) 

The main results are 



s-i ({-i)e + (o-i)A 2 3 , 

V2,Hopf ~ -JZTl^ 1 (Q _ 1)3 Pl + ^l)) 



5-1 90(6 - l) 3 + 5A(9 - l) 3 2 ^, 3 . 
= {e _ m25e _ 5 _ e) Mi + O^), 

lj_ = -C [9 (6(6 - 1)9 + 9(9 - 1)A)] , C> 0. 
For the Hopf bifurcation curve we impose the conditions x = 0, y = and §§ + §| 



0. Solving a series expansion yields the result for the curve. Next, the first Lyapunov 
coefficient l\ is computed using the invariant formula (5.39) from [23] that leads to (2.10). 

For the heteroclinic curve we proceed as follows [5] . We assume 5, 9 < and 59 — 1 > 
and we transform variables to obtain a system that is a perturbation of a Hamiltonian 
system. This enables us to formulate a Melnikov function. Setting this function to zero 
yields an equation from which we extract the quadratic approximation to the heteroclinic 
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curve. We introduce the transformation (t, x, y, ^2) — > {ex p ~ 1 y q ~ 1 t,ex,ey,—e,cie + 
C2e 2 ) where 

5-1 1-5 1-9 
Cl ~ I^T ,P ~ 56^T ,q ~ 59^1' 

Then we obtain 

V J \ Vipy + 5x-y) J y \ Ac 2 y + yx 2 J ' 

which for e = is a Hamilton system with Hamiltonian 

H(x, y) = -x p y q (-l + x + j^yj ■ 

Define gi = Qx p y q+1 and §2 = x p ~ 1 y q (02 + As 2 ) . The Melnikov function along H(x, y) = 
h is given by the following integral 

M(h) = / gidy-g 2 dx 

JH=h 

= I -Ox p y q+1 dy - x p - 1 y q (c 2 + Ax 2 ) dx 

JH=h 

= [ -(x p - l y q+2 ^- + x p ~ l y q (c 2 + Ax 2 ))dx, 

where we used Green's Theorem to convert the dy term to dx. Now along the nontrivial 
critical curve H(x, y) = we have y = — x ) so that 

"«■> - - (j^jjy^ ~ +A ) * 



Qp{5 - I) 2 
( q + 2)(9-lf 



C2 i P -i, q + e v li9+2 ( j- ;\ (Q \ + a ) , 



where we defined 



Jo r(2 + a + b) 



Solving M(0) = and substituting p, q we obtain 

« = v-mws-,) iee(s - 1,3 + 4A( " - 1)3 > • 

As a final check we consider the difference between the Hopf and heteroclinic curves 

(50-l)h 2 3, 

H2,HET ~ ^HOPF = ce{9 _ m2§9 _ § _ fl + 0(A*l). 

We see that the quadratic approximations of the curves coincide when the Hopf bifurcation 
is degenerate. 

Appendix C. Higher order coefficients 

In this appendix we list the third order normal form coefficients for the LPNS bifurcation 
and the fourth and fifth order coefficients for PDNS and NSNS, which are necessary to 
determine the stability of the tori (if they exist). Remark that we have not listed the 
coefficients (or necessary coefficients of the expansion of the critical center manifold) which 
can be obtained by complex conjugacy or the similar expressions for u>2 instead of uj\ in 
the case of NSNS. 
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C.l. Third order coefficients for LPNS. The normal form coefficients in (2.3): 
1 f T 

a-soo = - / , C(r; v\, v\, Vi)+3B(t; «i, /i20o)+3/i20o-6a 2 oo^200-6a200^(T)ui) aV+a 20 o 
6 Jo 

1 f T 

b 2 io = o / («2> C( r ; vi,v 1 ,v 2 )+B(t; v 2 , h 200 )+2B(T; v 1 ,h no )+2h lw -2a 2 ooA(T)v 2 ) dr+b 1 
1 Jo 

i r 

^021 = r / (v 2 , C(t; v 2 , v 2 , U2) + -B(r; u 2 , ^020) + 2.B(t; u 2 , /toil) - 2a iiA(r)i;2) dr 
^ Jo 

a m = / (</>*, C(t;v 1 ,v 2 , v 2 ) + 2$t(B(T;v 2 ,h 10 i)) + B(t;vx, h u) + /ton - aoa-AOr)^ 
Jo 

-23?(6no)/ioii - aon«20o) <^t + aon 

C.2. Fourth and fifth order coefficients for PDNS. The normal form coefficients in 
(2.4): 

1 r T . 

0400 = 24 y ^ )- D ( r ) - ' u i ! ' y i)' u i)' u i) + 6C f (r;ui,«i,/i 2 oo) + 3B(r; /t 20 o, ^200) 
+4B(r;«i,/i 3 oo) - 12a 20 o/t20o) = 



"211 = 77/ (V 7 *,- ^; ^1,^1,^2,^2) + C(r; ui,ui,/iou) + C(r;u 2 ,U2,/t20o) +4K(C(r;vi,v 2 ,/n 



+2K(B(r; u 2) /i 20 i)) + B(t; h 200 ,h on ) + 2B(t; h m ,h no ) + 2B(t; v 1: h n i] 
-"011^200 - 2q 2 oo«oii) dr 



1 



T 



"022 = j I (<p*,D(t;v 2 ,v 2 ,v 2 ,v 2 ) + 4C(T;v 2 ,v 2 ,h ii) + 2$l(C(T;v 2 ,v 2 ,hoo2)) 







+£(r; h 020 , hoo2) + 2B(t; hou, h u) + 45?( J B(r; u 2 , /i i 2 )) 
-4aonAon) dr 



Fourth order coefficients of the expansion of the critical center manifold can be computed 
by solving the following BVPs on [0,T]: 

A 400 - A{r)h m o ~ D(t;v 1 ,v 1 ,v 1 ,v 1 ) - 6C(r; vi, vt, h wo ) - 3B(t; h 200 , h 200 ) 

-4S(t;ui,/i 3 oo) + 12a 2 oo/t200 + 24a 4 oow + 24a 30 o/t200 = 

hioo(T) - /t 4 oo(0) = 

f T (<P*,h 4 oo) dr = 

/i 40 - M T )how + 4«w/i 40 - D(t; v 2 ,v 2 ,v 2 ,v 2 ) - 6C(r; v 2 ,v 2 , h 02 o) 

-AB(T;v 2 ,ho3o) - 3B(t; h 020 , h 02 o) = 

h 0i0 (T) - h 040 (0) = 
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^310 - ^(t)^310 + iuh 3 io - D(t; vi,vx,vi,v 2 ) - 3C(r; vi,vi,hn Q ) 
-3C(T;vi,v 2 ,h 2 oo) - B(T]V2 : h 300 ) - 3B(T;vi,h 2 ia) - SB(t; ^200,^110) 

+6a 2 oohno + 6a 30 o/ino + 66210^110 + 6ia;a 2 oo^no = 

h 3 io(T) + h 310 (0) = 

h 130 - A(r)h 130 + 3iujh 130 - D{t\ v 1 ,v 2 ,v 2 ,v 2 ) - 3C(r; v 2 ,v 2 , hno) 
-3C(T;vx,v 2 ,h 02 o) - B(T;vi,h 03 o) - 3B(t; h 020 , h 110 ) - 3B(T;v 2 ,h 120 ) = 

hisoiT) + /»i3o(0) = 

^031 - A(r)h 031 + 2iu;h 031 - D{t\ v 2 ,v 2 ,v 2 ,v 2 ) - 3C(r; v 2 ,v 2 , h 011 ) 
-3C(T;v 2 ,v 2 ,h 020 ) - B(T;v 2 ,h 030 ) - 3B(T;h 020 ,h 011 ) - 3B(T;v 2 ,h 021 ) 

+3a n/io20 + 66021/1020 + 6iwa n^020 = 
h 03 i(T) - h 031 (0) = 

h 2 u - A{r)h 2 \\ - D(t;v 1 ,vi,v 2 ,v 2 ) - C(T;v 1 ,v 1 ,h 011 ) - C(T;v 2 ,v 2 ,h 200 ) 
-m(C(T;v 1 ,v 2 ,h 10 i)) - 2$t(B(T;v 2 ,h 201 )) -£(r; /i 20 o, /tpn) 
-2B(t; hioi,hno) - 2B(t; vi,hm) + a 0U h 200 + 2a 2 oo^on 

+2a 2 nu + 2a m /i 2 oo + 43?(6 2 io)^on = 
/i 211 (T)-/i 211 (0) = 

Jo((p*,h2ii)dT = 

hi2i ~ A(r)h 12 i + iuh 121 - D(t; v 1 ,v 2 ,v 2 ,v 2 ) - C(r; vt,v 2 , h 02 o) 
-2C(r; vi,v 2 , h Q11 ) - C(r; v 2 , v 2 , h 101 ) - 2C(r; v 2 , v 2 , hno) - B(t; v 1 ,h 02 i) 
-B(T;h 020 ,h wl ) - 2B(T',hpu,huo) - 2B(r;v 2 ,hiu) - B(T;v 2 ,h 120 ) 

+2a n/ino + 26 02 i/ino + 2a m /ino + 2ia;a ii/ino = 

hi2i(T) + hm(0) = 

^220 - A( T ) h 220 + 2iuh 220 - D(t;v 1 ,v 1 ,v 2 ,v 2 ) - C(T;v 2 ,v 2 ,h 200 ) 
-4C(r; v 1 ,v 2 , h lw ) - C(r; «i, «i, /1020) - #(t; ^200, ^020) - 2S(r; v 2 , h 2W ) 
-25 (r; hxio, h 110 ) - 2B(t; vi, h 120 ) + 2a 200 h 020 + 46 2 i /i 20 + 4zwa 20 o/i020 = 

h 220 (T) - h 220 (0) = 

h 022 - A(T)h 022 - D(t; v 2 ,v 2 ,v 2 ,v 2 ) - 4C(r; v 2 ,v 2 , h n) - 25R(C(t; v 2 , v 2 , h 002 )) 

-B(T;h 020 ,h 002 ) - 2B(t; /ton, /ion) - 45ft(5(r; v 2 , h 012 )) 

+4a;oiiAoii + 4a 22-uo + 83?(6 2i)fy)ii = 
h 022 (T) - /io2 2 (0) = 

fo{f*,h 022 )dT = 

Fifth order normal form coefficients in (2.4): 



a 500 



1 /" T 

— y (u 1 ,£;(T;?;i,?;i,vi,'i;i,?;i) + 10D{t; «i, «i, «i, /i 2 oo) + 10C(r; vi, vi, /i 30 o) 

+15C(r; «i, /i 200 , ^200) + 10S(t; /i 20 o, ^300) + 55(r; «i, /i 4 oo) - 20a 20 o^300 
-120a 40 o-4(r)fi) dr - a 20 oa 30 o 
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foio = J {v2,E(t;vi,vi,v 1 ,vi,v 2 ) + 6D(T;vi,vi,v 2 ,h 2 oo) + 4:D(T;vi,vi,v 1 ,h 110 ) 

+4C(r;vi,f2,/i30o) + 6C(r; v±, v±, h 2 w) + 3C(r; v 2 , fo o, fo 00 ) + 12C(r; vi, h 200 , h lw ) 
+45(t; vi, h 3W ) + 4£(r; foio, fooo) + 6B(t; fooo, foio) + B(t; v 2 , fo 00 ) 
-24a i0Q A(T)v 2 - 12a 20 ofoio) dr - a 20 o&2io 



a 311 = - (v*1,E(t;v 1 ,vi,vi,v 2 ,v 2 ) + D(T;vi,vi,vi,h 011 ) + 6%l(D(T;vi,vi,v 2 ,h 101 )) 
o Jo 

+3D(t; V!,v 2 , v 2 , fooo) + 3C(r; fooo, foil) + 6K(C(r; u 2 , fooo, fooi)) 

+6K((7(t; vi,v 2 , h 201 )) + C(r; u 2 , v 2 , h 300 ) 

+3C(r; «i, «i, /im) + 6C(r; «i, fooi, foio) + 3£(r; fo 00 , fon) 

+6^(B(t; h 201 ,h 110 )) + 2K(B(t; v 2 , fooi)) + B(t; foil, fooo) 

+3B(r;/i 2 n,fi) - 6a 2 nA(r)wi 

-"onfooo - 6a 2 oofon) - 0200^111 - aona300 



&221 = jj (v 2 ,E(t;v 1 ,v 1 ,v 2 ,v 2 ,v 2 ) + D(T;v 2 ,v 2 ,v 2 ,h 200 ) 
+2D(T;v 1 ,v 2 ,v 2 ,h 101 ) + 2D(T;v 1 ,v 1 ,v 2 ,h 011 ) 

+D(T;v 1 ,v 1 ,v 2 ,h 020 ) + 4£)(r;i;i,'U2,«2,^iio) + 2C(r; u 2 , foio, foio) + C(r; ui, ui, fo 2 i) 
+C(r; u 2 , u 2 , fooi) + C(t; v 2 , fo 00 , h 020 ) + 2C(r; u 2 , u 2 , foio) 

+2C(T;vi,v 2 ,hi2o) + 2C(r;ui,/i 20,^ioi) + 4C(r; «i, v 2 , ftm) + 4C(r; t> 2 , fooi, foio) 

+2C(r; v 2 , fooo, h 011 ) + 4C(r; vi, h 110 , foil) + #(t; v 2 , fo 2 o) 

+2 J B(r; ui, fo 2i ) + 2B(t; fo 2 o, fooi) + AB(t; h 110 , foil) + 2B(t; h 2W , h 011 ) 

+2B(t; v 2 , fon) + B(t; fooo, fo2i) + B(t; fooi, fo2o) 

-2aonfoio - 4q 2 h^1(t)i; 2 - 2a 20 ofo2i) dr - 0200^021 - "011^210 



0122 = 4y (vI,E(t;vi,v 2 ,v 2 ,v 2 ,v 2 ) + 4^(D(T]v 2 ,v 2 ,v 2 ,hi i)) + 2^(D(T-,vi,v 2 ,v 2 ,h o 2 )) 

+4£>(r; vi,v 2 ,v 2 , foil) + 2C(r; vi, fon, fon) 

+8K(C(r; w 2 , fon, fooi)) + 4K(C(r; w 2 , fo 02 , foio)) 

+C(r; vi, fo 2 o, foo2) + 4K(C(r; ui, u 2 , foi2)) + 4C(r; u 2 , u 2 , fon) 

+2K(C(r; u 2 , u 2 , h 102 )) + 4S(t; fon, fon) 

+2K(£(t; fo 20 , fo 02 )) + 4SR(£(r; v 2 , fo i2 )) + B(t; «i, fo 22 ) 

+4K(B(r; foio, foi2)) - 4a nfon 

-4q 22-4(t)-ui) dr - aoiiam 



PERIODIC NORMAL FORMS FOR CODIM 2 BIFURCATIONS OF LIMIT CYCLES 



11 



fy)32 = To I ( v 2, E ( T ; v 2,V2,V2,V2,v 2 ) + D(t; v 2 , v 2 , v 2 , h 002 ) + 3D(t; v 2 , v 2 , v 2 , h 02 o) 
11 Jo 

+6D(T;v 2 ,v 2 ,v 2 ,h 011 ) + 6C(t; v 2 , h 20, ^011) + 6C(r; v 2 , v 2 , h Q21 ) + C(r; v 2 , v 2 , 7i 3o) 
+3C(r; v 2 , h 002 ,h 020 ) + 3C(r; u 2 , v 2 , h 012 ) + 6C(r; v 2 , /toil, ^Oli) 
+3B(T;h 020 ,h 02 i) + 6B(t; h 02 i, /toil) + 3-B(r; «2, ^022) + 2-B(r; u 2 , /to3i) 
+5(t; h 002 , h 030 ) - 6a nh 02 i - ^ao22M T ) v 2) dr - "011^021 

C.3. Fourth and fifth order coefficients for NSNS. Fourth order normal form coef- 
ficients for (2.5): 

1 f T 

"2200 = 4y ,%t)i,Di,!)i,vi) + 2SR(C(r;vi,«i,/iQ2oo)) 

+4C(r; «i, vi, ftiioo) + -B(t; /t 20 oo, ^0200) + 25(r; hnoo, ^1100) + 43?(.B(t; vi,hi 200 )) 
— 4anoo^iioo) dr 



ami = / {^*,D(T;vi,vi,V2,V2) + C(T;vi,vi,hooii) + 2St(C(T;vi,v 2 ,ho 10 i)) 
Jo 

+2^{C(T;v 1 ,v 2 ,h no)) 

+C(r; v 2 ,v 2 , hnoo) + 2K(5(r; «i, h iu)) + B(t; Ttoiio, ^1001) + B(r; h wi, h 1010 ) 
+B(t; /toon, hnoo) + 25ft(5(r; v 2 , h 1101 )) 
— "ooii^noo — anoo/toon) dr 



Fourth order coefficients of the expansion of the critical center manifold can be computed 
by solving the following BVPs on [0,T]: 



hiooo ~ A(T)h moo + 4iwi/i 40 oo - D(t;v 1 ,v 1 ,v 1 ,v 1 ) - 6C(r; v\, t>i, /t 20 oo) 

-3B(r;/t 20 oo,/t 2 ooo) - 4B(r;v l5 /13000) = 
^4000 (r) - /t4ooo(0) = 

/t3ioo - ^(r)/t3ioo + 2mji/i 3 ioo - £>(r; «i, «i, «i, «i) - 3C(r; t>i, v 1: hnoo) 
-3C(r; «i, wi, /t 2 ooo) - B(r; «i, /13000) - 31? (r; «i, /t 2 ioo) - 3,B(t; /i 20 oo, ^noo) 

+3anoo/t20oo + 6a 2 ioo/t20oo + ^101100^2000 = 

/*3iooCT)-/»3ioo(0) = 



how ~ A(T)h 30 io + 3iujih 3010 + iw 2 /i 30 io - D(t; vi,vi,vi,v 2 ) - 3C(r; V\,v\, /11010) 

< -3C(r;t;i,t; 2 ,/i 2 ooo) - -B(t;u 2 ,/i 30 oo) - 3-B(r;ui,/t 2 oio) - 3-B(r;/t 20 oo,/tioio) = 

^30io(T)-/t 3 oio(0) = 

ft300i -^4(T)/t 3 ooi +3^i/t 30 oi -W2/13001 - ■D(t;«i,«i,«i,«2) - 3C(t;«i,i;i,/hooi) 

< -3C(r;t;i,'U2,/t200o) - -B(t;u 2 ,/i 30 oo) - 3-B(r;ui,/t 2 ooi) - 3-B(r;/t 2 ooo,^iooi) = 

Tt300l(T)-/l 3 00l(0) = 
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h 22 oo - A(r)h 2 2ao ~ D(t; v 1} v 1 , v\, v\) - 2^(C(T;v 1 ,v 1 ,h 02 oo)) 
-4C(r; vi, vi, h no o) - B(t; /i 20 oo, ^0200) - 2B(t; h uoo , h noo ) - 4M(B(t; vi, /11200)) 

+83?(a 2 ioo)^noo + 4a 2 20o'"o + 4anoo^iioo 
h 22 oo(T) - /i 2 20o(0) 

A2020 - A(T)h 2 o20 + 2iwi/i 2 o20 + 2w 2 h 2020 - D(r;ui,«i,u 2 ,?; 2 ) - C(r; vi, vi, /10020) 
-C(r; w 2 , v 2 , h 2000 ) - 4C(t; vi,v 2 , frioio) - #(t; /i 20 oo, ^0020) - 2B(t; u 2 , fr 20 io) 

-2B(t; frioio, feioio) - 25(r; wi, /11020) 
7*2020 (T) - /l202o(0) 

/12002 -^4(t)/i 2 oo2 +2zwi/i 2002 - 2«w 2 /i 20 o 2 - -D(r; vi,vi,v 2 ,v 2 ) - C(r; u 2 , u 2 , ^2000) 
-4C(r;t>i,u 2 ,/iiooi) - C(t; «i, «i, ^0002) - 25 (t; v 2 , 7i 2 ooi) - B{r; h 20 oo, h 0002 ) 

-25(r;/iiooi,/iiooi) - 2£(t; ui, /11002) 
h 2 oo2(T) - h 2002 (0) 

h 2 no - A(r)h 211 Q + zwi/i 2 no + iu 2 h 2 u ~ d (t; vi,vi,v~i,v 2 ) 
-C(r; t>i, «i, hollo) - 2C(t; t>i, «i, /iioio) ~ G(t; «i, «2, ^2000) - 2C(r; «i, v 2 , ^1100) 
-S(r;«i,/i 2 oio) - 25(t;/iioio,/iiioo) - B(t;v 2 , /i^ioo) - #(r; ^-2000, ^0110) 
-2B(t; vi, hm ) + 2a 2 i o/iioio + 26ino/iioio + 2anoo^ioio + 2i(u>i + w 2 )anoo^ioio 

h 2 no(T) - /i2iio(0) 

^2101 - -4(t)/i2101 + wi/i 2 ioi - iuj 2 h 2W1 - D(t; vi,vi,vi,v 2 ) - C(r; v\, vi, h ioi) 
-2C(T;vi,vi,hiooi) - C(r; v\, v 2 , /t 2 ooo) - 2C(T;v 1 ,v 2 ,h noo ) - 2B(t; h 100 i, /inoo) 
-2B(r;«i,/inoi) - -B(r;v 2 ,/i2ioo) - B{t\v\, /12001) - -B(t; /t 2 0O0, ^0101) 
+2a 2 ioo^iooi + 261101^1001 + 2aiioo^iooi + 2i(a/i — w 2 )anoo^iooi 

^2101 (T) - /i2ioi(0) 

h 2 ou - A(T)h 2011 + 2iwi/i2oii - D(t;v 1 ,v 1 ,v 2 ,v 2 ) - C(r; vi, v%, /i on) 
-2C(r; ui, «2, ^1001) - C(t; t> 2 , ^2, ^2000) - 2C(r; ui, %, /11010) - 5(r; w 2 , /i 20 io) 
-B(r; v 2 ,/i 2 ooi) - B(t; h 200Q , h 00 n) - 2B(r; /tiooi, ftaoio) ~ 2#(r; v%, /lion) 

+2aion/i20oo + aooii^-2000 + 2iwiaoon^20oo = 

^oncn - ^2011(0) = 

Ann - A(r)h nil - D(r;v 1 ,v 1 ,v 2 ,v 2 ) - C(T;v 1 ,v 1 ,h oii) 
-23?(C(t; vt,v 2 , h 0101 )) - 2K(C(r; vi,v 2 , ^0110)) 
-C(r; w 2 , U2, ^1100) - 25R(B(r; ui, /10111)) - B(r; h 0110 , h 100 i) 
-B(t; h i 01 , h ww ) - B{t\ h Q011 , huoo) - 25R(B(r; v 2 , h U01 )) 

+2K(aom)/iiioo + "0011^1100 + 2K(6noi)/ioon + aim^o + anoo^-0011 = 

hnu(T) - h nn (0) = 

Jo (^^mi) dr = 

Fifth order normal form coefficients for (2.5): 
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■T 



O3200 = T7T / {vI,E(t;v 1 ,vi,v 1 ,vi,vi) + D(T;v 1 ,v 1 ,v 1 ,h 02 Qo) + 3D(T]Vi,V 1 ,Vi,h 2 OOo) 



1Z Jo 

+6D(t; vi,v 1 ,vi,h 1W0 ) + 6C(r; «i, /in 00 , h 1100 ) + 3C(r; ui, vi, ^1200) 
+C(T;v 1 ,v 1 ,h 30 oo) + 6C(r;wi,t;i,/i2ioo) + 6C(r;wi, fr 2 ooo, ^1100) + 3C(r;ui, /i 200, ^2000) 
+-B(r; /i 200, ^3000) + 2-B(r; vi, h 3W0 ) + 3£(r; ui, /12200) + 6B(r; h 2 ioo, ^1100) 
+3£>(r; /1200O) ^1200) — 6anoo/i2ioo 
-12a>220oA( T ) v i) dT ~ "110002100 



&0032 = 777/ (v 2 ,E(t;V 2 ,V 2 ,V 2 ,V 2 ,V 2 ) + f 2,^2,^2, /i0002) + 3-D(t; t> 2 , ^2 , ^2 , /lQ02o) 



^ Jo 

+6D(T;v 2 ,V2,V2,h 0011 ) + 6C(r; v 2 , h 0011 , h 0011 ) + 3C(t;v 2 ,v 2 , h 0012 ) 

+C(t;v 2 ,v 2 , /ioo3o) + 6C(t; u 2 , u 2 , ^0021) + 6C(r; -y 2 , /10020, ^0011) + 3C(r; w 2 , /iooo2, /10020) 

+B(t; /i oo2, /ioo3o) + 2B(r; v 2 , h 0031 ) + 3B(t; v 2 , h 0022 ) + 6B(t; h 002 i,h 0011 ) 

+3B(t; /ioo20, ^0012) — 6aoon/ioo2i 

-12a 022-4(T)w2) dr - aoon^xm 



01022 = J (vI,E(t;v 1 ,V2,V2,V2,v 2 ) + D(t;vi,v 2 , v 2 , ^0020) + D(t;v 1 ,v 2 , v 2 , /10002) 



+2£>(r; u 2 , u 2 , «2, ^1001) + 2£>(t; v 2 ,v 2 , v 2 , h ww ) 
+4D(T;v 1 ,v 2 ,v 2 ,h 0011 ) + 2C(t;vi,v 2 , h 0(m ) + C(r; ui, /i 00 20, ^0002) 

+2C(r; ui, u 2 , ^0012) + 2C(r; u 2 , ^1001, %)2o) + 4C(r; u 2 , /iiooi, ^0011) + 2C(r; u 2 , /iioio, ^0002) 
+C(r; v 2 ,v 2 , h W02 ) + 4C(r; u 2 , /iioio, ^0011) + 4C(r; u 2 , v 2 , h wu ) + C(r; u 2 , v 2 , h 1020 ) 
+2C(r;ui, /loon, ^0011) 

+ J B(r; ui, /ioo22) + 25 (r; ^0021, ^1001) + B(t; h 0020 , h W02 ) + 2B(t; h 0012 , h ww ) 
+4B(t; /loon, ^1011) + 2£(r; u 2 , fc.1012) + #(t; ^0002, ^1020) + 2-B(r; u 2 , ^1021) 
-4a oii^ion - 4a o22-4(7")-ui) dr - aooiiaioii 
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■T 



h2io = 7/ (v2,E(t;v 1 ,v 1 ,v 1 ,v 1 ,v 2 ) + D(T;v 1 ,v 1 ,v 2 ,h 2 ooo) + D(T;vi,v 1 ,v 2 ,h 0200 ) 



* Jo 

+2D(T;v 1 ,v 1 ,v 1 ,h 0110 ) + 2D(T;v 1 ,v 1 ,v 1 ,h ww ) 

+4D(t;v 1 ,v 1 ,v 2 , huoo) + 2C(T;v 1 ,v 2 ,h 2W0 ) + C(t;v 2 , h 2000l h 0200 ) 

+2C(T;v 1 ,v 2 ,h 1200 ) + 2C(t;v!, h 0110 , h 2000 ) + 4C(r; vi, h 01w , hi 100 ) + 2C(t;v 1i h ww ,h 020 o) 
+C(r; v 1 , V!, h 0210 ) + 4C(r; ui, h 1010 , h 1W0 ) + 4C(r; «i, vi, h llw ) + C(r; ui, v 1 , h 20W ) 
+2C(T;v 2 ,h 1100 ,h 1100 ) 

+B(t; v 2 , h 2200 ) + 2B(t; h 2W0 , h 01w ) + B(t; h 2000 ,h 0210 ) + 25(r; huoo, ^ioio) 
+45(r; Ziuoo, ^1110) + 2B{t\v 1 , h 12 io) + B(t; h 0200 , h 20W ) + 2S(t; ui, /i 2 no) 
-4anoo^ino - 4o!2200-4(t)w2) dr - 0440004440 



«24i4 = 0/ (^,^(T;fi,fi,-U4,w 2 ,W2) + £>(r;ui,t;i,U2,/ioioi) + D(T;v 1 ,v 1 ,v 1 ,h on) 



* Jo 

+D(T;v 1 ,v 1 ,v 2 ,h uo) + 2D(T;v 1 ,v 1 ,v 2 ,h ww ) + 2D(r; ui, v 2 , v 2 , h 1100 ) 
+D(t;v!,v 2 ,v 2 , h 2000 ) + 2D(T;vi,vi,v 2 ,h 1001 ) + 2C(t;v 1 , h 1001 , h 0110 ) 
+C(t; v 2 , v 2 , h 2100 ) + 2C(r; v 2 , h W01 , h 1100 ) + 2C(r; v ly h 1100 , h oou ) 
+2C(T;v 1 ,h wl o,h 0101 ) + C(r; «i,v 2 , /i20io) + C(r; vi, v 2 , /i 20 04) 
+2C(r;ui, v 2 ,/i mo ) + 2C(t;i1 2 , foioio, ^1100) + 2C(r; -ui, v 2 , h 1101 ) 
+C(r; u 2 , /i 2 ooo, h 0W i) + 2C(r; ui, /iioo4, h W i ) + C(r; ^1, /i 20 oo, fy)on) 
+C(r; u 2 , /i 20 oo, /i044o) + 2C(r; vi, vi, 04041) + C(r; V!,vi, h 0111 ) 
+B(t; v 2 , h 2W1 ) + B(t; h 2W0 ,h 0011 ) + 2B(t; vi, h lul ) + B(t; v u h 2011 ) 
+B(t; h 0101 , h 20W ) + B(t; h 2001 ,h 0110 ) + B(t; h 2000 , h 0111 ) + B(t; v 2 , h 2110 ) 
+2B(t; h wll ,h 1W0 ) + 2B(t; h ww , h 1101 ) + 2B(t; h 100 i,h 1110 ) 
-2a>nnA(T)vi - 2an o04044 - «oo44^240o) dr - 0004402400 - "440004044 
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&1121 = 0/ (v2,E(t;v 1 ,v 1 ,v 2 ,v 2 ,v 2 ) + D(t;v 1 ,v 2 ,v 2 , hoioi) + D(t;v 2 ,v 2 ,v 2 , hnoo) 
1 Jo 

+D(t; vi,v 2 ,v 2 , h WO i) + 2D(t; v 1 ,v 2 ,v 2 , h ww ) + 2D(t; v 1 ,v 1 ,v 2 , h 0011 ) 
+D(t; v 1 ,v 1 ,v 2 ,h o2o) + 2D(T;v 1 ,v 2 ,v 2 ,h 0110 ) + 2C(r; v 2 , h 0110 , h WO i) 
+C(T;v 1 ,v 1 ,h 0021 ) + 2C(r;t;i,/ioiio,^ooii) + 2C(r; v 2 , h 0011 , h 1100 ) 
+2C(t; v 2 , h ww ,h 0101 ) + C(r; vi,v 2 , h W20 ) + C(r; «i, v 2 , h 0120 ) 
+2C(r; «i, u 2 , /iion) + 2C(r; «i, /iioio, ^0011) + 2C(r; vi,v 2 , h 0111 ) 
+C(r; ui, /ioo20, ^0101) + 2C(r; v 2 , /10110, ^1010) + C(r; u 2 , ^0020, ^1100) 
+C(r;ui,/i 020,^iooi) + 2C(T;v 2 ,v 2 ,h 1110 ) + C(t;v 2 ,v 2 , h 1101 ) 
+B(t; V!,h 0121 ) + S(r;/i 02i,^iioo) + 2B(t;v 2 , h ull ) + B(r; u 2 , ^1120) 
+B(t; h 0W1 ,h W20 ) + B(t; h 0120 , h WO i) + B(t; h 0020 , h lwl ) + B(t; ui, ^1021) 
+25(r; /lino, /toon) + 25(r; /ii io, /loin) + 25(r; /lono, ^1011) 
-2a>nnA(T)v 2 - 2aoon/iino - "1100^0021) dr - an o^oo2i - "001161110 
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